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ANALYSIS OF ELECTROMAGNETIC INTERFERENCE FROM 
POWER SYSTEM PROCESSING AND TRANSMISSION 
COMPONENTS FOR SPACE STATION FREEDOM 


1 Introduction 

The goal of this research project is to analyze the potential effects of 
electromagnetic interference (EMI) originating from power system processing and 
transmission components for Space Station Freedom. The approach consists of four 
steps: 

1. Develop analytical tools (models and computer programs). 

2. Conduct parameterization (what if?) studies. 

3. Predict the global space station EMI environment. 

4. Provide a basis for modification of EMI standards. 

This report summarizes work performed during the period January 1 - August 
31, 1991. 

2 Device Models 

This portion of the investigation centers on the development of computer 
models which can be used to predict the local electromagnetic fields for various power 
system components or devices, including various transmission line configurations. 
Recent work is summarized in the Master’s thesis of Mr. Zemin Luo, which is 
attached as Appendix A. 

3 Radiation Sources 

This portion of the investigation centers on the development of computer 
models to simulate the propagation of electromagnetic waves in the ionospheric 
plasma and the radiation from antennas simulating EMI sources on the space station. 
Recent work is summarized in the Master’s thesis of Mr. Brian J. Hurysz, which is 
attached as Appendix B. 

4 Future Work 

Continuing work involves an investigation of vx B effects and the development 
of a more sophisticated method-of-moments model of the space station. 
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Abstract 


Ballooning techniques in conjunction with the two dimensional finite element (2DFE) 
method are described for modeling and computing a two dimensional (2D) magne- 
tostatic field in a two-wire parallel-pair transmission line problem in an unbounded 
region. Solutions using the combined ballooning 2DFE method are shown in this 
thesis to agree closely with analytical solutions for the two-wire parallel-pair trans- 
mission line in the same unbounded region. These 2DFE techniques are used here 
to study the effects of magnetic shielding and the effects of defects in such shield- 
ing. The distribution and profiles of magnetic fields surrounding such transmission 
lines are studied, for future use of investigating electromagnetic compatibility aspects 
associated with power processors and systems in the future Space Station Freedom. 
The combined ballooning 2DFE method is developed further for a 2D diffusion type 
(eddy-current) field and then applied to solve the problem in the two-wire parallel- 
pair transmission line with a conductive shield. Various cases such as the effects of 
conductive shielding thickness, and the effects of defective shielding are considered. 

Furthermore, in this thesis, a closed-form solution for the magnetostatic case in a 
two-wire twisted-pair transmision line is derived using Biot — Savart Law. The flux 
density comparisons between twisted-pair and parallel-pair cases are carried out and 
the effects of the length of the pitch of the twisted-pair line on the magnetic field 
profiles are obtained. The results shown in this thesis imply that a three-dimensional 
(3D) ballooning method using finite elements is needed for the magnetic shielding 
case as well as in the quasi-static diffusion type case. This quasi-static class can 
involve transmission lines in plasma, with or without conductive metallic shielding. 
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Chapter 1 

INTRODUCTION 


This work is one portion of a large research project centering on the “Analysis of 
Electromagnetic Interference from Power System Processing and Transmission Com- 
ponents for Space Station Freedom ”, and is concerned with aspects of developing 
device models for power system processing components and transmission lines in 
aerospace, where an infinite boundary and a plasma environment exist. 


1.1 Background and Literature Search 

One of the purposes of this work is to develop capabilities to compute the electromag- 
netic fields surrounding two-wire parallel-pair transmission lines and other possible 
transmission line configurations in aerospace applications, which can be solved as two 
dimensional (2D) problems with infinite (open) boundaries. 

The finite element method was applied to magnetostatic field solutions which 
involved the calculation of electromagnetic fields in bounded regions governed by 
quasi-Poissonian partial differential equations (PDE) in 1970 [1]. Boundary condi- 
tions involving Neumann, Dirichlet and periodicity types were involved [1]. While 
such boundary conditions occur when magnetic fields are confined by magnetic-type 
or conductive-type materials, there are many problems without such confinement. 
Many field problems are unbounded; i.e., they cannot be encased within a finite 
boundary along which either the normal or tangential component of the field may be 
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specified. As a consequence, the actual boundary is at infinity, with the potential zero 
thereupon. When dealing with integral methods, this boundary condition is automat- 
ically taken care of and a greater economy of nodes and elements are given, however, 
the integral methods have disadvantages in that the resulting matrices are dense and 
their formulation and inversion are much less straightforward than the equivalent dif- 
ferential approach. However, it is difficult to account for partial differential equation 
(p.d.e.) inhomogeneities with differential methods for an infinite boundary. On the 
other hand, differential methods can easily account for PDE inhomogeneities, while 
requiring that the boundary be nearby so that the global matrix for the bounded 
region solution is not excessively large. Approximating an infinite (open) boundary 
by a closer one makes the potential fall (or rise) more rapidly and distorts the solution 
from the true open boundary results. 

In a large number of electric and magnetic field problems, the solution requires 
the representation of a boundary which is remote from the region of interest. Typical 
of this class of problems are the single-sided linear machine, bus-bar impedance, 
and so on. There is, consequently, a need for a technique to combine a differential 
approach in the region containing all the field sources and materials with a simple 
representation of the exterior empty space. Thus additional work by Lowther (1978) 
[3], and Antunes (1982) [2], on a recursive (ballooning) finite element technique was 
developed for 2D planar magnetostatic field problems with open boundaries extending 
to infinity. Brauer (19S2) [4] developed that ballooning technique for axisymmetric 
problems and general skin effect problems in which the eddy current is confined in 
the interior region and the Laplacian equation is satisfied in the exterior region. 

Sometimes one is interested in the fields in both an interior region and an exterior 
region in a given problem. However, the above mentioned ballooning technique only 
uses the exterior region for the solution of the interior region. The solution of the 
exterior region cannot be obtained by the ballooning technique. Overcoming this dif- 
ficulty with exterior field problems has spurred a whole new area of research activity; 
as a result, several schemes are now available for tackling this aspect. Some types 
of algorithms of semi-analytic methods (or hybrid harmonic/finite element method) 
tackling open boundary problems were presented by Lee and Cendes (19S7) [5], Chari 
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(1987) [6], and Chari and Bedrosian (1987,1988) (7, 8] respectively. They used a cir- 
cle as a boundary. The solution in the exterior region is obtained analytically and 
matched to the interior finite element solution. 

To study the electromagnetic fields surrounding a two- wire parallel-pair transmis- 
sion line in aerospace applications, a ballooning technique can be used because we are 
only interested in the fields in the interior region near the wires. However, no such 
ballooning techniques have been found in any literature for that type of problems 
because eddy current fields can exist in both the interior region and exterior region in 
the expected plasma environment of Space Station Freedom. This thesis will derive 
a new 2DFE ballooning method for the problem with eddy currents in both the inte- 
rior region and exterior region and with an open boundary. Then this new ballooning 
2DFE method is applied in this research which will be defined further in the next 
section. 

The other purpose of this work is to develop models which can help engineers 
in studying the electromagnetic fields surrounding a two-wire twisted-pair transmis- 
sion line and possible other configurations. To the best of this author’s knowledge, 
no similar papers on such problem have been found in the literature. It is one of 
the purposes of this work to explain the nature of the field caused by this type of 
transmission lines. In this thesis, we only consider the magnetostatic case for the 
two-wire twisted-pair transmission line without shielding so that the closed-form an- 
alytical formulation can be employed to deal with this possible three dimensional 
(3D) problem. The results of this investigation will show the necessity of the need for 
the ballooning in conjunction with 3DFE methods to solve the eddy current problem 
with a twisted-pair transmission line with shielding. 


1.2 Definition of the Problem 

This research work centers on the development of computer models which can be used 
to predict the local electromagnetic fields with infinite boundary for various power 
system components or devices in aerospace applications, and specifically the case of 
the Space Station Freedom. 
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The effort is started in the magnetostatic case and is developed in the quasi- 
static case because of the existence of a plasma environment surrounding the Space 
Station. The first power system component to be analyzed is the two- wire parallel-pair 
transmission line with no shielding. This example is chosen because it is the simplest 
case of a two- wire transmission line for which there exists an analytical solution for 
comparison with and verification of the detailed FE modeling which will follow. 

A two-dimensional finite element (2DFE) method will be used to obtain the lo- 
cal electromagnetic field distributions surrounding power system components and 
devices. In order to simulate the problem with the open boundary, a ballooning tech- 
nique is introduced. With the ballooning technique, a matrix recursion algorithm is 
derived for a 2D problem in the static case and incorporated into the 2DFE method. 
This combined 2DFE ballooning method is applied to the example. The validity of 
this 2DFE method with the ballooning technique is demonstrated by comparison of 
the results with those of the analytical solution in the magnetostatic case with no 
shielding. 

In the magnetostatic case, the 2DFE ballooning method is applied to the two- 
w’ire parallel-pair transmission line with continuous permeable cylindrical shields to 
determine the effects of the shielding. This method is applied further to the study of 
the field distributions resulting from discontinuitues in permeable shielding surfaces 
due to accidental damage or to deliberate hardware design necessities. 

The 2DFE ballooning model is developed in this work for the AC cases with 
frequencies in the Hz.and kHz ranges. In the AC case, the magnetic vector potential 
is a function of time. Assuming that the MVP varies sinusoidally, we can use phasor 
vectors in the problems. The field in the plasma environment can be equivalent to 
the eddy current field. A complex 2DFE model for eddy current problems is obtained 
and a complex matrix recursion algorithm is derived using the ballooning technique. 
In this thesis, it is shown how the AC ballooning is working in each chosen frequency 
(/= 0 Hz, 400 Hz, 2 kHz, and 20 kHz). With the AC 2DFE ballooning method, the 
effects of conductive shielding are studied in lower frequencies and higher frequencies. 
The effects of defective conducting shielding are also considered in a certain frequency 
(/= 400 Hz). The relations between the flux density magnitude and shield thickness 
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are shown in both defective shields and shields without defect cases. 

The investigation of a 3D closed-form analytical formulation using the Biot — 
Savart Law is depicted for the magnetostatic field surrounding the two- wire twisted- 
pair transmission line. That formulation is an integral expression and the Gaussian 
numerical integral method is needed for the computer program. The results of this 
investigation indicate whether the near-field or the far-field are of 3D nature. Several 
different lengths of pitches are chosen for the flux density component comparisons, 
which also show that this problem cannot be simplified by a 2D solution approach. 
The 3D nature of this field means that it will be necessary to use 3DFE analysis to 
evaluate the magnetic material-type or eddy current-type electromagnetic shields and 
other related studies. 
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Chapter 2 

THEORETICAL 

BACKGROUND 


2.1 Formulation - the Governing Differential Equa- 
tions 

In this chapter the differential equation governing the magnetic fields within a two 
dimensional (2D) continuum is derived from Maxwell’s equations and the definition 
of the magnetic vector potential. For the purpose of this study, the magnetic fields 
surrounding the two- wire parallel-pair transmission line are treated as being two di- 
mensional (2D) in nature of the infinitely long transmission line case. Consequently, 
the conductor current can have a component in the axial-direction only. The materials 
within this continuum are assumed isotropic. 

The fields in such a device are governed by Maxwell’s equations, equations (2-1) 
through (2-4), as follows [11]: 


< 

tol 

II 

o 

(2.1) 

V-D = p 

(2.2) 

— - - dD 


VxH = J, + J e + — 

(2.3) 
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(2.4) 


Vxl = 


dB 

dt 


All the symbols used throughout this thesis are defined in the list of symbols (Ap- 
pendix A). The constitutive relationships are: 


B - fiH = -H (2.5) 

v 

D = tE ( 2 . 6 ) 


J t = uE 


(2.7) 


The magnetic vector potential^ and the magnetic flux density B are related as: 


B = V x A 


( 2 . 8 ) 


The x and y components of A are zero in the 2D case, therefore B can be rewritten 
as follows: 


-o _ dA *~ dA *- 

B ~ dy az dx° v 

Rewriting the curl of A in terms of the magnetic field intensity gives: 


(2,9) 


i/B = H = i/(V x A) 


( 2 . 10 ) 


Substituting Equation(2.10) into the third Maxwell equation, Equation(2.3), one ob- 
tains: 


V x H = V x (i/V x A) = J. + J e -f — 

oi 

In 2D cartesian coordinates Equation(2.11) becomes: 
d , dA t s , d , dA t , dD t 

aJ ( ‘ / a7 ) + = ~ {J - + J - + "aT* 


( 2 . 11 ) 


( 2 . 12 ) 


In the static case, the eddy current density component, J e , and the displacement 
current density component, dD/dt are both equal to zero, hence the contributions 
due to these terms will be zero. In the diffusion (eddy current) case, the displacement 
current is very small compared to the source current. For this reason, the contribution 


i 


of the displacement current can be neglected. Assuming that the magnetic vector 
potentials are functions of time which vary sinusoidally in both the diffusion (eddy 
current) case and high frequency radiation case, one can use phasor form equations 
in the frequency domain in these two cases. For convenience the subscripts of A z and 
J z will be dropped in future developments and discussions throughout this thesis. 
Equation(2.12) can therefore be rewritten as follows in the static case: 


d . dA d . dA 


(2.13) 


In the diffusion (eddy current) case, one can write the following: 


d , dA. d.dA = . ; 


(2.14) 


Also, in the high frequency radiation case, one can write the following: 


d . dA. d . dA. ? • i 2 ; 

+ = + (215) 

These are the well known nonlinear Poisson’s equations in cartesian coordinates. 

Usually, the reluctivity, 1 /, is a function of both position and magnetic flux density 
which is dependent upon the curl of the magnetic vector potential. However, in this 
w T ork, no nonlinear magnetic materials are expected to be encountered. Hence, v is 
assumed as a constant in the various anticipated nonmagnetic materials. Accordingly, 
the above three non-linear Poisson’s equations become linear equations. In other 
words, the reluctivity, y, can be moved outside of the partial derivative operators. 

Obviously a closed form solution is nearly impossible to solve the magnetic fields 
surrounding a two-wire parallel-pair transmission line with any shielding. One must 
look towards numerical techniques . One of the best methods is the finite element 
method which will be briefly summerized in the following section. 


2.2 The Finite Element Method 

The finite element method is used to solve bounded problems. Consider a 2D magnetic 
field problem with Dirichlet and Neumann boundary conditions. In the static case, 
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one can formulate the problem as follows: 

f ,d 2 A d 2 A, 

I ^ dx 2 + dy 2 ^ “ ~ Jt 


y4 = 0 

dA 

^ = ° 


( Governing p.d.e.) 
( Dirichlet ) 

{N eumann) 


( 2 . 16 ) 


From variational principles it was previously shown that the appropriate energy 
functional [ll] for this problem is: 

F(A) = UM ( fy?^ ( 2 . 17 , 

Minimization of F(A) is equivalent to solving the partial differential equation (p.d.e.) 
with Dirichlet or Neumann boundary conditions [11]. The solution region R can 
be divided into many subregions (finite elements) on which interpolation functions 
of the MVP’s are defined. There are many different types of element shapes and 
interpolation functions which are applicable to the finite element discretization [11], 
Here, first order triangular elements R e are chosen and interpolation functions can 
be obtained in the form of polynomials of first order. The value of the MVP at any 
point P(x,y) within this triangle is given by [11]: 

A e (x,y) = N,A, + N m A m + N n A n (2lg) 

where A’,,A’ m , and A'„ are the form (or shape) functions which are complete first 
degree polynomials. Meanwhile, A„A m , and .4. are the values of the MVPs at the 
vertices of the triangle, e, in Figure(2.1). These polynomials are of the form: 


( 2 . 19 ) 


the constants a,, b„ and c, are determined subject to the following conditions: 

N < = >. N m = N n = 0 at P(xi,y,) 

A " = >• A'i = A'. = 0 at P(x m ,y m ) 

— 1, A r m = A) = 0 at P(x n .y T) ) 


' N, ‘ 


ai 

bi 

Cl 


’ 1 ' 

N m 

= 

G-m 

bm 

Cm 

• 

X 

K 

^ J 



bn 

C„ 


. y . 


( 2 . 20 ) 
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Therefore, the expressions of the form (or shape) functions can be rewritten as follows: 
N, pi qi r, 1 

= 2^ Pm q m r m x (2.21) 

_ iV„ J p n q n r n J [ y _ 

where: A is the area of the triangular element, e, and 
n = In - Xm,qi = Vm - J in, Pi = X m y n - X n J/ m 
r m = Xj - x„, q m = y n - yt,p m = - x,y„ 

^ n = *£/ } === J// J/m^Pn “ 

where, 2 A = pi 4- p m + p n - 

Now, one can write the interpolation function of the MVP in one element, e, as 
follows: 

= 23 ( p > + Q'* + r *y) A ' ( 2 - 22 ) 

ts/,m,n 

Substituting y4 e (x,y) into the previous functional equation and minimizing it, one 
can get the following equation in matrix form [11]: 

S?A' = F (2.23) 


In expanded form we have: 



where: 


S-j = ^( r * r i + n 

T AJ * ■ 1 

/, = -T- t = /,m,n 

It should be noted that 
1 1 1 
•£/ ^"n 

y/ ym y n 



( 2 . 24 ) 

( 2 . 25 ) 

( 2 . 26 ) 

( 2 . 27 ) 
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The element matrices j£ f given in the above summary are only valid for element 
e. Representation of the entire region R requires that these element matrices be 
assembled into a global coefficient matrix £. Then a global matrix equation can be 
obtained as follows, after imposition of the proper boundary conditions: 

£-A = l (2.28) 

The assembly process can be briefly summarized as follows: 

Step(l): Initialize all entries of the £ and ]_ matrices to zero. 

Step(2): Add each contribution of an element matrix to its corresponding location 
in the global j£ matrix. 

Step(3): Add each entry of the element matrix £ e to its corresponding location in the 
global forcing function, i.. 

Step(4): Repeat step(2) and step(3) until all the elements have been assembled, 
now’ the system of Equation(2.28) would be ready for the numerical solutions. 
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Figure 2.1: The First Order Triangular Finite Element 
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Chapter 3 


COMPARISONS BETWEEN 
THE CLOSED-FORM 
SOLUTION AND 2DFE 
SOLUTION FOR THE 
TWO-WIRE PARALLEL-PAIR 
CASE 


In this chapter, the magnetostatic fields surrounding a two-wire parallel-pair trans- 
mission line which is shown in Figure(3.1) are computed by the analytical and 2DFE 
methods. The results are compared for verification purposes. 


3.1 The Closed-Form Solution 


According to Ampere's circuital law: 

in ■ dl = I 
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(3.1) 


one can get the following equations for conductor #1 in Figure(3.1): 

jH^dh = 1 

r2* 

/ H\r\da\ = / 

Jo 

Hence, 

Hi • 27rr x = I 
Therefore, 

i 

Hi = 


27TT i 


That is, 

Hi, x = -Hi sin(aj) = - 


Hi, v = Hi cos(qi) = 


l y 

27T?'i Ti 


H’ 


X — 


2*ri n 

Similarly, one can obtain the equations for conductor #2, Figure(3.1): 

I y 


H 2, x = H 2 sin(a2) = 


27TT2 T 2 


H2, y = — H 2 cos(q 2 ) = 


w 

I * + T 


2 nr 3 r 2 

With superposition of the fields due to the two conductors, the expressions 
total flux density can be obtained as follows: 

»1 ,\ 1 


B x = g-y(-j - -2) 


2i r 


W 


I x-— x + — 

B = — ( 2_ _ _ 

v 2 7 T ( r? 


W 
*-) 


(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 
of the 

(3.10) 

(3.11) 
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Because of the definition: 


V xA = B 


(3.12) 


in the 2D case, equation(3.12) yields the following: 


B x = 


dA 

dy 


By = 


dA 

dx 


(3.13) 

(3.14) 


Substituting equation(3.10) into (3.13) and integrating yields the following result: 

.4 = “ ln(— ) (3.15) 

27t t i 

the quantities rj and r 2 are defined in Figure(3.1). With the above formulations, the 
MVP and flux density were computed for the case of Figure(3.1). The equipotentials 
near to the conductors are shown in Figure(3.2) and the equipotentials away from the 
conductors are shown in Figure(3.3). The MVP values and the flux density values at 
various locations were computed and are listed in Tables(3.1) through (3.4). 


3.2 The 2DFE Solution 

Let us consider the region shown in Figure(3.4), which has a natural open boundary 
near to the conductors. In conductor # 1, a current of 1=500 A flows in the direction 
out of the plane of the paper, and in conductor # 2, a same value of the current 
flows in the opposite direction. Its FE grid is shown in Figure(3.5). This bounded 
problem can be solved using the 2DFE method summarized in the previous chapter. 
The equipotentials obtained from the FE method with the grid of Figure(3.5) are 
shown in Figure(3.6). One also needs to investigate the distribution of the magnetic 
fields away from the conductors. Taking a region which has a remote open natural 
boundary as shown in Figur(3-7), we can solve this bounded problem using the 2DFE 
method and obtain the equipotentials shown in Figur(3.8). In the case of Figure(3.7), 
some MVP values and flux density values at various locations are listed in Tables(3.1) 
through (3.4). 
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3.3 Comparisons Between the Closed-Form Solu- 
tion and 2D-FE Solution 

In the two above sections, the analytical method and the 2DFE method were used 
to solve the same problem. The “true” solution obtained from the analytical method 
was compared with the solution obtained from the 2DFE method. If the solution 
obtained from the 2DFE method is not within reasonable bounds in comparison to 
the true closed form analytical one, means of improvement in the 2DFE method, and 
solution algorithm need to be arrived at. 

Comparison between Figure(3.2) and Figure(3.6) leads one to conclude that the 
equipotential lines obtained from the 2DFE method are obviously different from the 
closed form solution. And comparing between Figure(3.3) and Figure(3.8), one finds 
that in the area near to the conductors, the equipotential lines obtained from the 
analytical method and 2DFE method are not much different, but in the area close to 
the boundary, the lines are heavily distorted. The value comparisons for the MVP’s 
and flux densities are given in Tables(3.1) through (3.4). From these tables, it is seen 
that the near-fields (close to the conductors) obtained from the 2DFE method have 
higher accuracy and the far-fields (close to the boundary) obtained from the 2DFE 
method have lower accuracy. In the other words, the more remote the boundary, the 
less difference the near-fields have from the analytical solutions. 

The facts from the above comparisons leads one to conclude that the farther 
away the boundary in a problem the more improvement are obtained in the near- 
field solution. In the enginerring problems, it is not possible and not necessary to 
include an infinite boundary. Accordingly, it is better to have a recursion algorithm 
to extend the boundary to as far away as possible, and to stop the recursion once no 
further additional improvement can be achieved. That is the essence of the ballooning 
technique which will be discussed in the next chapter. 
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Figure 3.1: The Two- Wire Parallel- Pair Transmission Line Geometry 
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Figure 3.3: Equipotential Contours Beyond the Vicinity of the Two- Wire Conductors 
Obtained From the Analytical Method 


Dimensions in inches 
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Figure 3.4: The Geometry of the Region on the Immidiate Vicinity of the Two-Wire Conductors 
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Figure 3.5: The FE Grid in the Immidiate Vicinity of the Two- Wire Conductors 
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Conductors Obtained From the 2DFE Method 
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Figure 3.7: The FE Grid Beyond the Immidiate Vicinity of the Two- Wire Conductors 
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Figure 3.8: Equipotential Contours Be 
Obtained From the 2DFE 
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Table 3.1: Comparison Between MVP’s Obtained from Analytical and 2DFE Methods 
at Locations Away From the Conductors 
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Table 3.2: Comparison Between Flux Densities Obtained from Analytical and 2DFE 
Methods at Locations Away From the Conductors 
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Table 3.3: Comparison Between MVP’s Obtained from Analytical and 2DFE Methods 
at Locations Close to the Conductors 
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Flux Densil 

ty (Wb/m 2 ) 

Difference 

% 
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Table 3.4: Comparison Between Flux Densities Obtained from Analytical and 2DFE 
Methods at Locations Close to the Conductors 





Chapter 4 

THE BALLOONING MODEL IN 
THE 2D MAGNETOSTATIC 
CASE 


Ballooning techniques to simulate infinite boundaries in 2D magnetostatic field so- 
lutions using finite elements were developed by earlier investigators in various appli- 
cations [11]. As seen in chapter 3, it was demonstrated that open boundary 2DFE 
solution of the field surrounding the two-wire parallel-pair transmisssion line in open 
space was inadequate. Hence, in this chapter ballooning techniques for 2DFE codes 
of the type used by this investigator will be detailed here. That is, in this chapter, a 
recursion algorithm will be derived using ballooning techniques [2, 3], and then this 
algorithm will be incorporated into the 2DFE method which was summarized earlier 
in chapter 2. This 2DFE ballooning model will be applied to an example and its 
results will be examined by comparison to the analytical solutions. With the 2DFE 
ballooning model, the study for the effects of the magnetic-type shielding and the ef- 
fects of defects in such shielding surrounding parallel-pair transmission lines is carried 
out. 
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4.1 Ballooning Technique in the Magnetostatic 
Case 


Consider the two round-wire parallel-pair transmission line shown earlier in Fig- 
ure(3.4), a finite element grid of which is shown here in Figure(4.1). Using con- 
ventional finite element techniques, the interior region in this FE grid which has a 
boundary To, may be divided into triangular finite elements as shown in the figure. 
Hence the node distribution on the surface To of the region R, Figure(4.1), is well de- 
fined. The FE global matrix equation can be obtained using the conventional 2DFE 
method for the region, and the boundary, To, as follows: 


Srr 

Srt 0 


Ar 


' Jr' 

_ Sr 0 R 

& 

o 

o 


Aq 


. 7 o . 


(4.1) 


where the subscript in denotes the interior region, the subscript To (or 0) denotes the 
nodes on boundary, To, and the subscript R denotes the nodes in the region R. 

Assume that the region exterior to R contains only free space, hence the external 
field is purely Laplacian. The first step in the representation of the exterior region is 
to define a finite element mesh for an annular region surrounding the region R. This 
annulus has certain properties which can best be summarized as follows: 

(a) The number of nodes on the inner and outer boundaries are the same and lie on 
lines radiating from a star point. 

(b) The outer boundary nodes for each successive ballooning annulus have a one-to- 
one correspondence with those nodes on the original solution region boundary, To. 
This leads to useful geometric similarity properties to the FE triangles in the bal- 
looning region. 

(c) There is a fixed ratio k between the radii of the outer and inner sets of nodes of 
every annulus, measured from the star point. 

Such an exterior annular region with boundary annuli To, Tj, ... is shown in Fig- 
ure(4.2). The exterior region should be matched with the interior region such as 
shown in Figure(4.3). 
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Hence, if the points on the inner boundary are denoted: 

A'< 0) = (x„,!/„) (0) 


(4.2) 


then the outer boundary nodes are given by 


*<" = (x„,y«) (1> = tXl 


( 0 ) 


(4.3) 


where k is the mapping ratio for the ballooning method at hand. 

The FE matrix, 5.„ . for this exterior annular region can now be constructed 
and a process of matrix condensation is used to remove any equations corresponding 
internal nodes occurring in the recursive ballooning step. 

The FE matrix equation for nodes on To and Tj becomes: 


= 


’Si 1 

S 12 


1 

O 

1 


O 

1 

. $21 

$22 _ 


1 

> 

1 


O 

1 


(4.4) 


where the subscript ext denotes the nodes on the exterior region, the subscript of 
MVPs, 0, denotes the inner nodes on the To and the subscript 1 denotes the nodes 
on the Ti. 

The contributions of each element on this annulus to are computed as follows: 


S-j = + Mj) m > n 


(4.5) 


where: A is the area of the triangular element, e, and 
n = En T'm = ^ n = &1 

Ql = J/m J/n i Qm = J in Vi i 2/m 

In the same way, a second annulus for the nodes on Ti and T 2 can be constructed: 


Z'A'* = 


’n l 

T'n 


‘ A x ' 


‘ 0 ‘ 

.n, 

Ti,. 


. A 2 . 


0 _ 


(4.6) 


Because this second annulus is geometrically similar to the first and has the same 
mapping ratio discussed earlier, it follows that one can write: 


X™ = kxi" = k 2 A'<°> 


(4.7) 


Applying the geometric similarity given in equation(4.7) yields a set of elemental 
coefficients, S£, equation(4.5), which remains unchanged as one progresses from one 
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ballooning annulus to the next. That means the matrix X 1 in equation(4.6) is equal 
to the matrix 5>^ . t in equation(4.4). That is, £ Tt is identical to X 1 , or 


z 1 = £x< 

The two annuli maj 7 now be combined to yield the following: 

0 


Tl 2 


(4.8) 


Sli 

5j 2 

*^21 

si 2 + n i 

0 

n 1 



1 — 

0 

■ 


1 

O 



II 

0 


— ■ 
1 


1 

0 

« 


(4.9) 


where A\ are nodal MVPs on the common boundary, Tj, between the two annuli. 
The equations for Ao, Ai, and Ai may be partially solved to eleminate A\, using the 
fact that 


A 1 = -(Sl 2 + Ti 1 )-'(S l 2l Ao + Tl 2 A 2 ) 

Substituting Equation(4.10) into Equation(4.9) leads to the following: 

s;, - s.ysj, + T!,)-'SI, -S' U (S' 71 + 




0,2 _ 


-TMS'n + T’,)- 1 ^ T] 2 - T],{S\ 2 a TA)- 1 ^ J 


hence, 
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xt 
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— T 2j ( S 22 + TA ) 1 •S'jj 

5A Sf 2 

L *$21 *^22 


(4.10) 


>io 

Ai 


(4.11) 


-SWS^ + T’.J-T,', 

Til ~ Tl,(S'n + J 


(4.12) 


The matrix equation (4.11) has the form of equation(4.4) but relates the two combined 
annuli, which now has a mapping ratio of fc 2 . 

From the (i-l)th recursive ballooning step, applying the method given in equa- 
tion(4.12) recursively, one could obtain the following equation: 
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which relates the MVPs, Ao and A t on the boundaries To and I\ of the ballooned 
solution region. 

For the (i)th application, one can obtain: 


1 = 


T« 

-Ml 

n%' 


Ai 


’ 0 ' 

T' 
L i 21 

m 
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(4.14) 


which relates the MVPs of the nodes on boundaries I\ and IYn . One must recall 
that the above mentioned geometric similarity property of the successive ballooning 
layers allow one to write the following: 


r 



(4.15) 


Combining equation(4.13) and equation(4.14), and eliminating A,, one can write 
the following main recursive relationship for the ballooning technique, which is used 
here: 
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That is, one can write the following for the ballooned exterior region: 
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(4.16) 


(4.17) 


which relates the nodal MVPs on the To and I\+i, and has a mapping ratio of k^ ,+i \ 
A new annular region can now be attached to the present outer nodes and the 
process can be repeated indefinitely. Hence a recursive system has been defined with 
the outer boundary moving away with the following geometric progression: 

(4.18) 


The system is rapidly convergent and the Su terms in equation(4.17) (representing 
the error in the exterior solution) quickly approaches zero. After m applications, one 
can write: 
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Finally, the S u terms are attached to the outer nodes of the interior solution region, 
i.e., the matrix equation(4.1) is combined with the matrix equation(4.19), yielding 
the following global FE system of equations: 

Srr 5/ir 0 0 A I 

Sr 0 R £r 0 r 0 + Su Sjj Aq ~ Io (4-20) 

0 SJ] S? 2 J [ A m J 0 _ 

where Aq are nodal MVPs on the boundary, To, and A m are the nodal MVPs on the 
boundary, r m . If r m is far enough, A m = A\r m — 1 * 0 and 5JJ — 1 ► 0. 

If A m = 0, and = 0, the previous equation(4.20) can be rewriten as follows: 

Srr Srt 0 A _ I 

_ Sr 0 R Sr 0 r 0 + SJ 7 ] _ A 0 _ Jo 

In the computer algorithm of the 2DFE ballooning method, m, the number of bal- 
looning layers (or annuli, which from now on will be used interchangably throughout 
this thesis), is chosen large enough to lead to negligible coefficients for the matrix, 
SJj. Specific values for m will be discussed in the next section, accompanied by some 
numerical results. 

The process of the 2DFE ballooning method can be summarized as follows: 
Step(l): Compute the global j£, n matrix of the interior region using the conventional 
2DFE method summerized in chapter 2, that is, obtain the FE global matrix £,■„ and 
the global forcing function / shown in equation(4.1). 

Step(2): Choose a proper m, the number of ballooning layers, and compute the con- 
tribution of the exterior annular region To and Tj to obtain the FE matrix £^ xi . 
Step(3): i=l 
Step(4): Let T! = 

Step(5): Use the main recursive equation(4.16) to combine the matrices £^ xt with £' 
and therefore, obtain the matrix . 

Step(6): If i is less than m, i = * + 1, go to step(4). 

Step(7): Incorporate the matrix 5^ into the FE global matrix £ in as shown in equa- 
tion^. 21). 

Step(8): Solve the FE equation(4.2l). 

Step(9): End. 


( 4 . 21 ) 
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4.2 Exploration of the 2DFE Ballooning Model 
by Test Examples 

Now this 2DFE ballooning model will be applied to an example of a two-wire parallel- 
pair transmission line problem with an infinite boundary , which was shown earlier in 
Figure(3.1), and its validity should be examined by the comparison of its results with 
those of the analytical solution whose formulations were given in equations (3.13) 
through (3.15). 

The interior region is shown in Figure(4.4) where the current flowing in the right 
conductor has a value of 500 amperes and a direction out of the plane of the paper, and 
the same value of current flows through the left conductor in the opposite direction, 
that is into the plane of the paper. Its FE mesh is shown in Figure(4.5), where the 
mapping ratio is fc=1.3 and the start piont is at P{xo, yo) = (0,0). 

The equipotentials of the theoretical “true” solution obtained from the analytical 
method are shown in Figure(4.6). Figures(4.7) through (4.12) show the equipotentials 
obtained from the 2DFE ballooning method with ballooning layers (annuli), m=0, 2, 
5, 10, 15, and 20, respectively. The quantitive comparisons of the MVP’s and flux 
densities at arbitrarily selected locations shown in Figure(4.13) were carried out. The 
MVP comparisons between the analytical and 2DFE ballooning solutions in various 
locations for the ballooning layers, m=0, 2, 5, 10, 15 and 20 are listed in Tables(4.1) 
and (4.2). Also, the flux density comparisons between the solutions obtained from 
these two methods are also listed in Tables(4.3) through (4.5) 

Comparing Figure(4.6), which shows the analytical solution, with Figures(4.7) 
through (4.12), one can easily find that while the number of ballooning layers, m, is 
lower (m=0, 2, and 5, which show various FE ballooning solutions, see Figures(4.7) 
through (4.9), respectively), the equipotentials obtained from the 2DFE ballooning 
solutions are obviously distorted, meanwhile when m is increased (m=10, 15, and 
20, see Figures(4.10) through (4.12), respectively), the equipotentials of the 2DFE 
ballooning solutions are closer and closer to those obtained from the “true” solution. 
Again, from Tables(4.1) and (4.2), one can find that the errors of the MVPs were 
reduced while the number of ballooning layers was going up, for example, the errors 
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of the MVPs are in the range of (20 ~ 70)% with m=0, and in the range of (10 ~ 40)% 
with m= 5. For m > 10, the errors are below 10 %. In Tables(4.3) through (4.5), the 
errors of flux densities bear a similar relationship to the number of ballooning layers, 
m, to that relationship in the case of the MVP data explained above. Tables(4.1) 
through (4.5) show one that while the number of ballooning layers, m, is big enough 
(here, m > 10), the errors of both MVPs and flux densities are smaller and smaller 
up to a number of ballooning layers, at which such errors are not further improved 
by further addition of more ballooning layers. Such lack of continued improvement is 
due to the inherent numerical errors associated with the FE discretization within the 
FE grid of the internal region. That is, further improvements can only be obtained by 
finer FE mesh discretizations in the interior region. Therefore, one can say that this 
2DFE ballooning model is useful for magnetostatic problems with infinite boundaries. 


4.3 The Effects of the shielding and the Effects 
of Defects 

This 2DFE ballooning model can be used to investigate the effects of the material-type 
magnetic shielding on the electromagnetic field distributions in a two-wire parallel- 
pair transmission line case. 

Let the cylindrical shell shown in Figure(4.4) surrounding the two-wire parallel- 
pair transmission line be made of material with relative permeability fi T > 1.0. Finite 
element magnetic field analysis was applied using the ballooning technique to simulate 
the infinite free space surrounding the conductors (infinite boundary), with assumed 
cylindrical shell relative permeabilities of /< r = 10, 100, 1000, and 5000. This was done 
to study the effects of shielding on the magnetic field outside the shield. The equipo- 
tentials are shown in Figures(4.14) through (4.17), for fi r = 10, 100, 1000, and 5000, 
respectively. These should be compared with the no shielding case in Figure(4.12) 
for a qualitative appreciation of the effects of the shield’s relative permeability, p r , on 
the field distribution. A quantitive comparison of magnetic flux densities at arbitrar- 
ily selected locations outside the shield, shown in Figure(4.18), is given in Table(4.6). 
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Notice the substantia] effect of the shield when fi r > 100, the flux densities outside the 
shielding are very small because of the magnetic shielding with higher permeabilities. 

Also this 2DFE ballooning model can be used to study field distributions resulting 
from discontinuities in shielding surfaces due to accidental damage or to deliberate 
hardware design necessities, such as openings for electrical leads. Examples given here 
assume three possible cylindrical shielding shell discontinuities around the two-wire 
parallel-pair transmission line under consideration. These defects and their corre- 
sponding flux plots for defect angles, S = 10°, 20°, and 30°, are showm in Figures(4.19) 
through (4.24), respectively. The effects of the shielding defects were monitored at 
the arbitrarily chosen locations given in Figure(4.25). The magnetic flux densities are 
given in Table(4.7). The numerical results show how defects can adversely affect the 
intended benifits of such shields, which is manifested here in a substantial increase in 
the flux densities in element # 1305. 

Next the combined FE-ballooning technique is developed and applied to the case 
of shielding involving diffusion type reduced eddy currents in metalic shields for AC 
type transmission line currents. 


37 




3S 


Figure 4.1. Interior FE Grid With Boundary r 0 for the Parallel-Pair Transmission Line 
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Figure 4.2: The Nature of FE Grids in Bounded Exterior Regions for Ballooning Puroses 
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(») 2D-FE Solution Region without Ballooning (b) 2D-FE Solution Region with Mooning 
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Figure 4.5: FE Grid for the Two-Wire Parallel-Pair Transmission Line Case 
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Figure 4.6: MVP Equipotentials for the Parallel-Pair Case Obtained From the Analytical Solution 
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Figure 4.8: MVP Equi potentials Obtained From the FE Solutions With 2 Ballooning Layer: 
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Table 4.1: MVP Comparison Between Analytical and 2D-FE Ballooning Solutions for m=0,2 and 5 
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Table 4.2: MVP Comparison Between Analytical and 2D-FE Ballooning Solutions for m=10,20 and 30 
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Table 4.3: Flux Density Comparison Between Analytical and 2D-FE Ballooning Solutions for m=0 and 2 
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Solution (Approximation of Infinite Boundary) 
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Table 4.4: Flux Density Comparison Between Analytical and 2D-FE Ballooning Solutions for m=5 and 10 
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Table 4.5: Flux Density Comparison Between Analytical and 2D-FE Ballooning Solutions for m=20 and 30 
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Figure 4.13: Finite Elements Chosen for MVP and Flux Density Data Compari 
in Tables (4.1) Through (4.5) 
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Figure 4.14: The Parallel-Pair MVP equipotentials Obtained From FE Solutions With 
Ballooning (m=15) Including Shielding (/ t r = 10) 
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Figure 4.15: The Parallel-Pair MVPEcjuipotentials Obtained From FE Solutions With 
Ballooning (m=15) Including Shielding (/< r = 100) 
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Figure 4.18: Finite Elements at Which the Effects of Shielding on Flux Densities are Given in Table(4.6) i 
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Figure 4.19: Geometry of a Defective Shield Case (6 = 10°) 
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Figure 4.20: The Parallel- Pair MVP 

(m=15) Including a Defe 
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Figure 4.23: Geometry of a Defective Shield Case (<5 = 30°) 



P Equipotentials Obtained With Ballooning 
elective Shielding Case (/t r = 1000,6 = 30°) 
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Figure 4.25: Finite Elements at Which the Effects of Defective Shielding on Flux Densities are Given in Table (4.7) 
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Chapter 5 


THE BALLOONING MODEL IN 
THE 2D DIFFUSION 
QUASI-STATIC CASE 


The 2DFE ballooning model derived in the previous chapter can be used to solve the 
magnetic field problems with open boundaries only in the static case. For the space 
station applications of interest here, many field problems with infinite boundaries 
should be dealt with in the quasi static case, so that an AC 2DFE ballooning model 
is needed for the AC eddy current problems with open boundaries. In this chapter, the 
previous ballooning 2DFE model will be developed further in the frequency domain for 
use w r ith the quasi static cases at hand. The new model will be applied to the study of 
electromagnetic field distributions surrounding a two-wire parallel-pair transmission 
line with conductive shielding in the (Hz to kHz) frequency range. The effects of 
shielding thickness and the effects of defects are also investigated. 


5.1 The 2DFE Ballooning Model in the Quasi 
Static Case 

5.1.1 The 2DFE Model in the Quasi Static Case 
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From chapter 2, one should use Equation(2.14) as the governing differential equa- 
tion in the 2D quasi static case in the frequency domain. Phasor vectors and phasor 
matrices should be used. The corresponding MVP boundary value problem can be 
formulated as follows [11]: 

O' 1 A d 2 A 

v {-Q -2 + ~Q~ y i) = -J* + jvaA ( the governing p.d.e.) 

.4 = 0 ( Dirichlet condition on boundary C\) (5-1) 


dA 

l/ !h =0 


(Neumann condition on boundary C 2 ) 

The energy functional for the above problem is given by the following [11]: 

F(A) - / + )!) - 1 ' A + v^ dii y < 5 - 2 > 

Minimizing the functional in equation(5.2) is equivalent to solving the partial differ- 
ential equation (p.d.e.) with Dirichlet and/or Neumann boundary conditions, equa- 
tion^. 1). Using the FE method summarized in chapter 2, equation(5.2) can be dis- 
cretized and minimized. Consequently, a complex matrix equation can be obtained 
for every triangular element as follows [11]: 

!'A e = (£ e +J£ e )A e = f (5-3) 

where: 


Z = 


VS 

K, 

K 


Vf m 

Km 


K 

Kn 

Kn 


Vij = + q,q 3 ) 


i,j as /, m,n 


Hi = 


■ 

WL 

K'n ' 

A 

‘ 2 

1 

1 ' 

Wml 

W'mm 


= JU<T n 

1 

2 

1 

. W* nl 

W'nm 


1 

1 

2 


(5.4) 


(5.5) 


(5.6) 



' h ' 


’ 1 ' 



_ A7 S 


— 

An 

3 

1 


. I* . 


1 


i = /, m, n 


(5.7) 
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n = %n ^T715 
^ m = X/ In, 

^ n = *^m 

and, 

1 1 1 

X/ x m x n 

Vl Vm J In 



Ql — 2/m J/n 
9m = Vn - Vl 
qn = 2// " S/m 


(5.8) 


(5.9) 


Assembling the element matrices £ and forcing function vectors ]_ given in Equa- 
tion^. 3) for all elements, one can obtain an FE complex global matrix equation as 
follows: 


RA = l 


(5.10) 


where: 

ne ne 

£ = E^' = D£' + h:') (si') 

C = 1 C = 1 

where symbol ne denotes the total number of elements. 

5.1.2 The Ballooning Technique in the Quasi Static Case 

In [4], Brauer presented a ballooning model to compute eddy current problems 
with open boundaries. In that model, eddy currents and all sources should be included 
in the interior region, hence the exterior region must be free space. In this thesis, a 
new AC 2DFE ballooning model will be derived for the diffusion type case in which 
the eddy currents can exist everywhere. In this new model, the total region is also 
divided into two subregions: the interior region which includes eddy currents and 
all sources ,and the exterior region which includes eddy currents only. The exterior 
region still consists of many annuli which have the same properties as in the previous 
chapter such as the start point, the fixed mapping ratio, and the same number of 
nodes on every boundary and the geometry similarity property to the FE triangular 
elements in the ballooning region of the magnetostatic case of chapter 4. 

This 2DFE ballooning model for the quasi static case is very similar to that for 
the static case discussed in the previous chapter. However, the difference between 
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these two models should be highlighted as given next. In the ballooning model for 
the static case, every annular matrix is the same, that is, X' = jSLo which leads to 
the simple recursive algorithm discussed earlier. In the quasi static case, the element 
matrices have imaginary components which are proportional to the element area, 
A. Therefore, one should compute the annular matrix X for every annular region 
bounded by I\ and r,+j. This annular FE contribution can be written as follows: 

FA"*' = 0 ( 5 . 12 ) 


where: 

na 

£' = £(£'+&"> < 5 - 13 > 

e=l 

where the symbol na denotes the number of elements in the annular region. The real 
components of the matrix X* are unchanged from one annular region to the next. 
Hence, one can obtain the following: 

j /‘ = £ ,_1 = £'~ 2 = . . . = £’ ( 5 . 14 ) 


However, the imaginary components of the matrix 


are not constant, that is: 



. A, 

3 “ a u 


' 2 

1 

1 ' 

. A jfc 2 ' 

’ 2 

1 

l 

1 

2 

1 


1 

2 

1 

1 

1 

2 _ 

1 

1 

2 


( 5 . 15 ) 


where, k is the mapping ratio between any two adjacent annuli. 

Now, the FE matrix, ^. x i , for the exterior region can be constructed and any 
internal nodal MVPs can be removed by a similar process of matrix condensation as 
shown in the last chapter. Combining the AC ballooning model and the AC 2DFE 
model, one can obtain the process of this AC 2DFE ballooning model as follows: 
Step(l): Using Equations(5.3) through (5.6) and (5.11), compute the complex FE 
global matrix and the forcing function 2 in the interior region /?, Figure(4.5). 
Step(2): Choose a proper m, the number of ballooning layers, and using equa- 
tions^.^) through (5.15), compute the matrix S c xt which is the contribution of 
the exterior annular region To and Tj. 

Step(3): i=l. 
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Step(4): Using Equations(5.12) through (5.15), compute the matrix £ which is the 
contribution of the exterior annular region T; and r,+j. 

Step(5) Using the relationship between the matrices and £ shown in equa- 
tion^. 16), obtain the matrix where: 



C<+1 Ci+l 
‘-’ll Oj 2 

C>+1 C'+l 
0 2 1 ‘-’22 


(5.16) 


Step(6): Repeat step(4) and step(5) until i > m. 

Step(7): Incorporate the ballooning matrix 5J 1 ) into the FE global equation (5.3), 
that is: 


Srr Srt 0 


A 


7 

_ Sro/? Sr 0 r 0 + ^Ti 


. ^ r ° 


„ A'o . 


Step(8): Solve for the complex MVPs from the equation(5.17). 
Step(9): End. 


5.2 The 2DFE Solution Using Ballooning Tech- 
niques with Conductive Shielding 

The AC 2DFE ballooning model for the eddy current case has been obtained in 
the last section. Now it can be applied to a two-wire parallel-pair transmission line 
problem with a conductive shield and an infinite boundary. The geometry of the 
interior region was shown in Figure(5.1), where a current which varies sinusoidally 
with a peak value 1=500 A flows through the right conductor and in the direction 
out of the plane of the paper, and the same value of current flows through the left 
conductor in the opposite direction, that is, into the plane of the paper. Eddy currents 
exist in the copper shielding and elsewhere is free space. For the ballooning model, 
here the start point is chosen at (xo,!/o) = (0,0), and k = 1.3 is chosen as the mapping 
ratio. The FE grid of the interior region is the same as shown in Figure(4.5). 

The computations were carried out using the AC 2DFE ballooning model derived 
in the previous section for the eddy current case in various chosen frequencies (/=0 
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Hz, 400 Hz, 2 kHz, and 20 kHz). Figures(5.2) through (5.6) show the equipotentials 
in the case of frequency /= 0 Hz, that is the magnetostatic case, and with ballooning 
layers, m=0, 2, 5, 10, and 20, respectively. Figures(5.7) through (5.11) show those 
case with m=0, 2, 5, 10, and 20, when /=400 Hz, respectively. Figures(5.12) through 
(5.16) show those in the case of /= 2 kHz and with ballooning layers, m=0, 2, 5, 
10, and 20, respectively. Figures(5.17) through (5.21) show r those in the case of 
/= 20 kHz and with the ballooning layers, m=0, 2, 5, 10, and 20, respectively. The 
magnitudes of flux densities at some selected locations shown in Figure(5.22) are 
plotted in Figures(5.23) through (5.25) for various conductive shielding thicknesses. 

Comparing Figures(5.2) through (5.21), one can arrive at some conclusions as 
follows: 

(1) In a certain frequency (here, in any one of the frequencies, /= 0 Hz, 400 Hz, 2 
kHz, and 20 kHz), while the number of ballooning layers, m, is larger and larger, the 
variation in the contours of the equipotentials is smaller and smaller, which means 
that when m is large enough, the equipotentials may almost converge to those of the 
“true” solution. 

(2) In lower frequency cases, such as /= 0 Hz, and /=400 Hz, whenever m is greater 
than 10, the equipotentials are almost unchangeable. In higher frequency cases, less 
ballooning layers are sufficient to obtain more accuracy of the equipotential contours. 
For example, in the case of /= 2 kHz, it was found that a number of ballooning layers, 
m= 5, is quite sufficient. It is interesting to point out that in the case of /= 20 kHz, no 
ballooning layers are needed. These results reveal that the higher the frequency, the 
strong the weakening effects of eddy currents in the shielding on fields outside that 
shielding. That is, at higher frequencies the fields outside the shielding are effectively 
diminished. 

(3) In Figures(5.23) through (5.25), three shielding thicknesses are used. Those figures 
show the effects of varying the conductive shielding thickness, from which one can 
find that the thicker the conductive shielding the more effective it is reducing the 
outer fields. That is, the weaker the fields outside the shielding. 
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5.3 The 2DFE Solution Using Ballooning Tech- 
niques with Defects of the Conductive Shield- 
ing 

In this section, the AC 2DFE ballooning model is used to study the field distribu- 
tions resulting from various defects in the conductive shielding. For shielding defect 
angles, 8 — 10°, 20°, and 30°, the corresponding FE solution regions are shown in 
Figures(4.19), (4.21), and (4.23), respectively. In the case of a frequency, /=400 
Hz, the shielding thickness, A 4 = 0.025 in, the number of ballooning layers, m=10, 
and an AC current with a peak value, 1=500 A is assumed to pass through the two 
conductors in the directions shown in the figures. The corresponding equipotential 
plots for angular defects, 8 = 10°, 20°, and 30°, are showm in Figures(5.26) through 
(5.28), respectively. The effects of the defects were monitored at the arbitrarily cho- 
sen locations given in Figure(5.29). In the case of a frequency, /=400 Hz, an angular 
defect, 8 = 20°, and a shield thickness, A 4 = 0.025 in, the magnitudes of magnetic 
flux densities at the locations shown in Figure(5.29) are given in Figure(5.30). 

In the case of a frequency, /=400 Hz, comparing Figures(5.26) through (5.28) 
which show the equipotentials for defect angles, 8 = 10°, 20°, and 30°, respectively, 
with Figure(5.11) which shows those with no defects in the shielding, leads one to 
conclude that the larger the defect angle, 8, the greater and greater the magnetic flux 
densities leaking in the neighbourhood of the defect area. Comparing Figure(5.30) 
which shows the effects of defects in the case of an angular defect, 8 = 20°, and a 
shielding thickness, A, = 0.025m, with Figures(5.23) through (5.25), which show the 
effects of the shielding thickness without deflects, one can see the slight degradation 
of the effectiveness of the magnetic field shielding near the defect region. 

The AC 2DFE ballooning model can be used to solve parallel-pair transmission 
line prblems with infinite boundaries. However, in space station applications, there 
are potential situations which might involve twisted-pair type transmission lines, as 
well as other configurations not available yet, for which the above model is not suit- 
able. One of the simplest methods to solve the simplest twisted-pair transmission line 
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problem is the closed-form analytical method using the Biot — Savart Law in con- 
junction with 3D integral type method. These will be discussed in the next chapter. 
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Figure 5.2: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions for 

f 0Hz With Conductive Shielding Without Ballooning - 
The Natural Open Boundary 
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Figure 5.7: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions for 
f = 400 Hz With Conductive Shielding Without Ballooning - 
The Natural Open Boundary 
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Figure 5.11: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions 1 
With Ballooning (m=20) for / = 400// 2 With Conductive Shielding 
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Figure 5.12: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions 
for / = 2 kHz With Conductive Shielding Without Ballooning - 
The Natural Open Boundary 
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Figure 5.16: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions 
With Ballooning (m=20) for / = 2 kHz With Conductive Shielding 



Figure 5.17: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions for 
/ = 20 kHz With Conductive Shielding Without Ballooning - 
The Natural OpenBoundary 
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Figure 5.20: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions 
With Ballooning (m=10) for / = 20 kHz With Conductive Shielding, 
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Figure 5.22: Finite Elements at Which the Effects of Shielding Thickness on Flux 
Densities are Given in Figures(5.23) Through (5.25) 
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Figure 5.23: Flux Density Magnitudes With Ballooning (m=10) and no Defects With 

a Conductive Shield Thickness, A s =0.1 in, for f=0 Hz, 400 Hz, 2 kHz, 

20 kHz and 100 kHz, for a Current 1=500 A 
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Figure 5.24: Flux Density Magnitudes With Ballooning (m=10) and no Defects With 
a Conductive Shield Thickness, A 3 =0. 05 in, for f=0 Hz, 400 Hz, 2 kHz, 
20 kHz and 100 kHz, for a Current 1=500 A 
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Figure 5.26: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions With Ballooning (m=10) 

for / = 400 Hz - the Case of a Conductive Shield With an Angular Defect, 6 = 10° Shield 
Thickness, A,,=0. 025 in, for a Current 1=500 A 
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Figure 5.27: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions With Ballooning (m=10) 

for / = 400 Hz - the Case of a Conductive Shield With an Angular Defect, 6 = 20° Shield 
Thickness,A„=0.025 in, for a Current 1=500 A 
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Figure 5.28: The Parallel-Pair MVP Equipotentials Obtained From FE Solutions With Ballooning (m=10) 

for / = 400 Hz - the Case of a Conductive Shield With an Angular Defect, 6 = 30°, Shield 
Thickness, A 5 =0. 025 in, for a Current 1=500 A 
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Densities are Given in Figure(5.30) 


Magnitude of Flux Density (Lines/in**2) 



Figure 5.30: Flux Density Magnitudes With Ballooning (m=10) and a Conductive 
Shield With an Angular Defect, 6 = 20°, Shield Thickness, A^O.025 in 
for f=0 Hz, 400 Hz, 2 kHz, 20 kHz and 100 kHz, for a Current 1=500 A 


Chapter 6 

CLOSED FORM SOLUTION OF 
THE TWISTED-PAIR 

I 

TRANSMISSION LINE 


In chapter 4 and chapter 5, the 2DFE ballooning models have been developed in 
the DC case and AC case and applied to the two-wire parallel-pair transmission line 
problems with infinite boundaries which might become one of the common problems 
in the space station power transmission line applications. Another potential problem 
is that of the magnetic field surrounding twisted-pair transmission lines, again with 
infinite boundaries. This will be shown here not to be amenable to solution by the 
2DFE ballooning models. Therefore, some new methods are needed for those twisted- 
pair problems. In this thesis, only the simplest twisted-pair case is considered so that 
closed-form formulations can be derived by the Biot — Savart Law coupled with 
integral methods. From the results of these closed-form solutions, some useful insight 
into the nature of the magnetic field surrounding such transmission lines will be 
obtained. 


6.1 The Derivation of the Formulation 

6.1.1 Biot — Savart Law 
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The magnetostatic field may be described by the magnetic vector potential (MVP) 
A , which, in a homogeneous uniform linear medium, is governed by [11]: 

V 2 A = -fiJ (6.1) 

Equation(6.1) is called the MVP’s Poisson equation and may be solved by the well 
known integral formulation as follows: 

J = JL[1± (6.2) 

47 t J r 

where, the symbol r denotes the distance from the field (or observation) point to the 
source point. The relationship between the MVP and the flux density is given by 

VxA = B (6.3) 


Premultiplying both side of equation(6.2) by the operator Vx, and using the defini- 
tion in equation(6.3), one could obtain an expression for the flux density as follows: 



Consider the following useful vector identity: 

V x (-) = - V x 7-7 x V(-) (6.5) 

r r v 

where, J is only a function of the position of the source point, which leads to Vx J = 0, 
1 7 

and V(-) = — Equation(6.4) can therefore be rewritten as follows: 
r r 3 

~B = — [ im-dv ( 6 . 6 ) 

47 - J r 

that is known as the Biot — Savart Law. 

If the cross section of the conductor is very small, the following relationship could 
be obtained: 

Jdv = Jdsdl = Idl (6.7) 


Using equation(6.7), one can rewrite equation(6.2) and equation(6.6) respectively, as 
follow's: 


A = 


fj.1 r dl 
4tt J r 


( 6 . 8 ) 
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and 


B = tL [ £*1 (6.9) 

An J r 

6.1.2 Analytical solutions for a two-wire twisted-pair transmission line by the 
Biot — Savart Law 

Now, consider a two-wire twisted-pair transmission line which consists of many 
pitches in free space shown in Figure(6.1). The details of the geometry in one pitch 
are shown in Figure(6.2), where a pair of DC currents (1=500 A) flow through the 
pair of spiralling conductors, that is, one is spinning along the positive z axis and the 
other along the negative 2 axis. 

Now, consider the MVP formulation. Using the integral MVP method, equa- 
tion^. 8), one can compute the MVP at any arbitary point, P(x,y,z), shown in 
Figure(6.2), as follows: 


A( x v g) = HSL f *1 

1 An J g r. An J* r 2 


dh 

9 r l 


fid 


dh 

r 2 


t!2i N y f A_il 

h 0 r l r 2 


(6.10) 


where: the symbol NP denotes the number of pitches, the subscript, s, denotes the 
line source, r 5 denotes the space vector from the source point on conductor #1 to 
the field (observation) point, r 2 denotes the space vector from the source point on 
conductor #2 to the field (observation) point. 

From the geometry shown in Figure(6.2), some useful relationships can be written 
as follows: 


f ii = r 0 cos 6 \ 


I/! = r 0 sin 0i 



{ 0j : (- nNP,nNP ) 


( 6 . 11 ) 


x 2 = r o cos0 2 = —ro cos 0i 
2/2 = r o sin0 2 = — rosin 0i 

* h ft 

22 = - w) = 

02 — 0\ “I" 7T 


( 6 . 12 ) 
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(6.13) 


r a = [(x - r o cos0i) 2 + [y - r 0 sin6>i) 2 + (z - ^i) 2 ]* 


t ~2 = [(x + r 0 cos <?i) 2 + {y + r o sin0i) 2 + (z - — 0i) 2 ]^ 


2x 

h_ 

2n ' 


dl x = — t 0 sin 6\d9\d x + r 0 cos 9\d9\d y + —d9\d z 

Z7T 


dl 2 = r 0 sin 6\d6\d x — r 0 cos + — d9\d z 

in 


(6.14) 

(6.15) 

(6.16) 


Substituting equations(6.11) through (6.16) into equation(6.10), one may obtain 
the MVP as follows: 

r o sin0ic?0i r o sin0id0i 


j = poi r 2n{i+l) 

: \TD/o '' 2 ‘XX 


4tt 


i=-NP/2 ' 


ro cos $\d9i r 0 cos % dO, 


r 2 


5fL_S?l 


(6.17) 


r i r 2 

Now, consider the flux density formulation. Using the Biot — Savart Law , equa- 
tion^. 9), one can compute the flux density at any arbitary point, P(x,y,z), shown 
in Figure(6.2), as follows: 


B 


VjJ f l dji x ^q/ r 1 $2 x 

4?r Jo rj 47 t Jo r\ 

fioI_ NP ^~ i [ l '+' dl\ X r° _ (fl 2 xr« 

4 * ,=irp /2 ^ 

Hol /■'•+! d?i x Fl dl 2 x r 2 


) 


_ _ ) (6-18) 

47r i=-NP/2 Jti r l r 2 

where and are the unit space vectors from the source points, Pi and P 2 , on 
conductor #1 and conductor #2 to the field (observation) point, P, which are given 
as follows: 

r° = — [(x - r 0 cos 9i)a x + {y - r o s\n0i)d y + (z - — Oi)a z ] = — (6.19) 

2?r r i 

= — [(x + r 0 cos Oi)a x + {y - r o sin0i)a y -f (z - —61 )a%] = — (6.20) 

r 2 r 2 
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Using vector operation rules, one may obtain the following: 


dli x r i = 


Qx Q>y &Z 

— r 0 sin 81 d 9 \ r 0 cos 6\ dOi — d0i 

2 *h 

x — r 0 cos y — ro sin 0i z — — 0 X 

It 


r 0 cos0id0i(z — —01 ) — {y — r 0 sin 0i)— dd x 

h h 

(x - r o cos0i)— d9i + r Q sin 9xd6i(z - — 0 X ) 

1 7T Z7T 

— r 0 sin 6id9i(y — ro sin 0i) — ro cos 6\d9i{x — ro cos 0 X ) 


( 6 . 21 ) 


and, 


-d'l 2 x r 2 = | -rosin 9 x d0x r 0 cos9id9 1 

x + r 0 cos 9 } y + r 0 sin z — — 0i 

h r 


r 0 cos9id9i(z — 77-^1) + (y + r 0 sin 0i)— d0i 

— (x 4- r 0 cos 0i )— ddi + r o sin0id0i(z — — 0i) 

It It 

— r 0 sin 8id0i(y + r 0 sin0i) — r 0 cos 9\d9\(x + r 0 cos 0i) 


( 6 . 22 ) 


Substituting equation(6.21) and equation(6.22) into equation(6.18), one can obtain 
the formulation of the flux density as follows: 


B-* 0 - 


NP /2- 1 ,2*(.+l) 


= t2i <2 f 

ri 


i=-NP / 2 ' 


r 0 cos 8i(z - ^-Bi) -{y- r 0 sin 01)77- 
(x - r o cos0i)— + r 0 sin0 x (z - — 0 X ) 
-r o sin0i(y — r o sin0i) — r 0 cos 0 x (x — r o cos0i) 



r o cos0i(z - — 0i ) + (y + r o sin0i)— 

ll T h ^ 7 T 

~(x + r o cos0i)— + r o sin0i(z - — 0i) 

Z7T Z7T 


1\ 


d0i (6.23) 


— r o sin0i(y + r 0 sin0 x ) — r o cos0i(x + r 0 cos0i) 


/ 
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6.2 The Computation of Flux Densities 


The closed-form analytical formulations of the magnetic field surrounding the two- 
wire twisted-pair transmission line, which is depicted schematically in Figure(6.2), 
was derived in the previous section. Those formulations can be programed with the 
numerical integral method (see Appendix B) for use on computers. The magnetic flux 
density, B, specifically the B x and B y components, can be computed at any cross- 
section and at any angle, 0, shown in Figure(6.2). Here, in this thesis, this calculation 
was performed at locations along the three concentric circles shown by the dotted lines 
in Figure(6.3). The resulting B field at various points along the three dotted circles 
is shown by the B vector arrow displays, whose direction indicates the orientation of 
the B field in the x — y plane or the z direction, and whose length is proportional to 
the magnitude of the resultant B field along the x — y plane or along the z axis. 

These B field results are shown in Figure(6.4) for the parallel-wire case discussed 
in chapter 4, and in Figure(6.5) and (6.6) for the twisted-pair case for a pitch, h = 5", 
and a plane with 6 = 0°. There are no 2 -component B fields resulting from the parallel- 
pair case as shown in Figure(6.4), while there are 2 -component B fields resulting from 
the twisted-pair case as shown in Figure(6.6). These B field calculations were repeated 
for a plane with 6 = 90°, and the results are shown in Figures(6.7) and (6.8). 

The effects of the choice of the pitch, h, showrn in Figure(6.2), on the flux den- 
sity components, B x , B v , B z and resultant magnitudes, ||£?||, are documented in 
Figures(6.9) through (6.12) at the mid-plane point for a transmission line of length 
20000 inches. For all these cases, in Figure(6.2), r 0 is equal to 1.0 inch, and the DC 
current flowing in each conductor is equal to 500 A. 

Figures(6.5) through (6.8) demonstrate the 3D nature of the B field surrounding 
a two-wire twisted-pair transmission line. Figures(6.9) through (6.12) show that 
shorter lengths of the twist pitch make the B fields outside the transmission line 
decay more quickly in the radially outward direction. That is, at locations far away 
from the transmission line, the 2 -component of the fields, B t , is of the same order of 
magnitudes of the x, and y components, and B y . Accordingly, the 2 -component, 
B z , could not be neglected in comparison to the B x and B y . That is, the twisted-pair 
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problem should be dealt with using 3D magnetic field computation methods. 
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Figure 6.3: Twisted-Pair Cross Section at 0 
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Figure 6.4: Flux Density Components B x and B y for the Parallel-Pair in the x — y Plane (I = 500.4 ) 

( No z-Component of Flux Density in This Case) 
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Figure 6.5: Flux Density Components D T and B y for the Twisted-Pair for a Pitch,/* 



120 


Figure 6.6: Flux Density Component B z for the Twisted-Pair for a Pitch, 
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• Figure 6.7: Flux Density Components B x and B y for the Twisted-Pair for a Pitch , h 
5 (=12.7 cm) in the 0 = 90° Plane (I = 500/1) 
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Figure 6.8: Flux Density Component B z for the Twisted-Pair for a Pitch, 
5”(=12.7 cm) in the 0 = 90° Plane (/ = 50(M) 



Plane: Z=0. L=20000 (in.) Theta=45 (deg.) 



Figure 6.9: Magnitude of Flux Density Component \B x \ max for Different Pitch 

Lengths (/ = 50(L4) 


123 




Plane: Z=0. L=20000 (in.) Theta=90 (deg.) 



Radius (in.) 


Figure 6.10: Magnitude of Flux Density Component \B v \ max for Different Pitch 
Lengths (I = 500.4) 
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Figure 6.11: Magnitude of Flux Density Component \B t \ m ax for Different Pitch 
Lengths (7 = 500.4) 
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Plane: Z=0. L=20000 (in.) Theta=90 (deg.) 



Radius (in.) 

Figure 6.12: Magnitude of Flux Density |B| max for Different Pitch Lengths 

(/ = 500 . 4 ) 






Chapter 7 


CONCLUSIONS AND FUTURE 
WORK 

7.1 Conclusions 

In this thesis, the ballooning methods incorporated into the 2DFE model have been 
presented to solve the 2D electromagnetic problems with infinite boundaries in the 
DC case and the quasi static (AC) case. The verification of this model was carried out 
by comparison of the results with those of the analytical solution in the magnetostatic 
case for a two-wire parallel-pair transmission line problem with no shielding and with 
an infinite boundary. This is because there exists a closed-form solution for that type 
of problem. Studies of the field distributions resulting from the magnetic shielding in 
the DC case and the conductive shielding in the AC case were performed, and the 
effects of possible defects in such shields were considered in both cases. 

In the application to solve two- wire parallel-pair transmission line problems with 
infinite boundaries, an interior region, I?, including eddy currents and all sources, 
with a circular boundary and an exterior annular region only representing free space 
were chosen. Usually, in this 2DFE ballooning model, not only the free space but the 
eddy current can be included in the exterior region, and because the field solution 
is limited to the interior region, /?, the consideration in choosing the contour To 
is that all regions of interest should be included in the interior region, R. This 
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2DFE ballooning model is applicable to complicated as well as simple shapes of the 
boundary, To, which must satisfy the criterion that no radial line from the start point 
Po intersects such a boundary more than once. 

Recalling the procedure of the AC 2DFE ballooning method, after modeling the 
interior region, R , by the FE method, it is necessary to form only one additional FE 
matrix, for an annulus one element wide which shares interface nodes with R, and 
then the FE matrix for the any other annulus is just the linear combination of the real 
components and complex components of the first one. With this performance of the 
annular matrices, the recursive computation to form the FE matrix for the exterior 
region becomes quicker and simpler. Having developed the recursive algorithm and 
obtained the FE matrix for the total exterior region, the only terms of interest in 
the matrix are those associated with the interface nodes on To, so that the size of 
the matrix which has to be solved to compute the field in the interior region, R, is 
unaffected by the addition of the exterior region. 

According to ballooning theory, the accuracy of the solution can be improved by 
increasing the number of ballooning layers, m. Recalling the results shown in chapter 
4 and chapter 5, one may find that the lower the number of ballooning layers, m, the 
more improvement can be achieved by increasing m. However, if m is already large, 
little improvement could be obtained by increasing its value. Actually, there are three 
facts which would affect the accuracy of the solution: (1) the number of ballooning 
layers, m; (2) the density of finite elements both in the interior region and the annuli; 
(3) the order of the interpolation polynomials (shape functions) used in the finite 
element formulation. Therefore, when the 2DFE ballooning model is applied to the 
magnetic field problems with infinite boundaries, the number of ballooning layers, m, 
should be chosen judiciously, and if further improvements of the results are required, 
finer FE mesh discretizations such as smaller finite elements, or high order finite 
elements will have to be adopted. 

A closed-form solution for the magnetic fields surrounding a two- wire twisted-pair 
transmission line in the static case was presented in the previous chapter. The results 
show that the shorter the pitch of the twist, the weaker the magnetic fields surrounding 
that twisted-pair transmission line. Comparing these results, one may find that three 


128 


components of the flux densities are in the same order of magnitude, hence, for near 
field effects none of them can be neglected. In other words, the twisted-pair problem 
is a 3D problem. 

7.2 Future Work 

According to the conclusions in the previous section, the 3D nature of the field sur- 
rounding a twisted-pair transmission line means that it will be necessary to use 3DFE 
analysis or some equivalent method to evaluate the effects of magnetic material-type 
shielding or eddy current-type electromagnetic field shielding in association with such 
twisted pairs. It may be possible and necessary to develop 3D ballooning models for 
twisted-pair problems or some other space station problems with infinite boundaries, 
which is one of the possible future works stemming from this research. 

In the space station applications, a future space station can be modeled as a long 
wire structure (with a diameter of 1 m and a length of 100 m approximately). In 
some cases, the space station orbits the earth at an altitude of 300 km with a velocity 
of approximately 8 km/s. The motion of the space station through the geomagnetic 
field results in an electromotive force v x B along its length. This electromotive force 
(emf) causes a potential drop across the space station and a current existing between 
the space station and the space plasma. 

A 3DFE method and a 3D ballooning method may be adopted to compute the 
electromagnetic field distribution for this problem. This may be accomplished as 
follows: 

(1) Choose a sphere as an interior region enclosing the space station. 

(2) Use a conventional 3DFE method to build up the global matrix system on the 
interior region where the v x B components are considered. 

(3) Divide the exterior region outside the interior region into many layers which are 
concentric with the interior region. 

(4) Compute the FE matrix for each layer and use ballooning techniques to obtain 
the contribution matrix for the exterior region. 

(5) Incorporate the contribution matrix of the exterior region into the global system 
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including the interior region to solve the problem. 

It seems possible to obtain the local electromagnetic field distribution in detail 
using such a 3DFE ballooning model. These fields can be used to calculate the 
currents, as well as physical quantities and phenomena of interest. 
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APPENDIX A 
List of Symbols 


Ay, A z 

A 

B X ,B V ,B Z 

B 

D 

E 

F(A) 

B 

Jxt Jyi Jz 

Js 

Je 

I 

V 

P 

a 

P 

t 

UJ 

s 

I 

A 

MVP 

p.d.e. 


The x, y and z components of the magnetic vector potential 

Magnetic vector potential 

The x, y and z components of the flux density 

Magnetic flux density 

Electric flux density 

The electric field indensity 

Energy functional 

The x, y and" z components of the field intensity 
Magnetic field intensity 

The x, y and z components of the source current density 

The source current density 

The eddy current density 

Source current 

Material reluctivity 

Material permeability 

Conductivity of medium 

Charge density 

Permitivity of medium 

Angular speed 

Global FE coefficient matrix 

Forcing function vector 

Vector of MVP components 

Magnetic vector potential 

Partial differential equation 
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AC Alternate current 

DC Direct current 

2DFE Two dimensional finite element 

3DFE Three dimensional finite element 
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APPENDIX B 

Gaussian Quadrature Integral Method 


There is a large family of numerical integration formulas that conform to the 
following pattern: 


J ^ f(x)dx » Aif(xi) + A 2 f(x 2 ) + . . . + A n f(x n ) 


In using such a formula, it is only necessary to know the “nodes” xi,X 2 ,...,x n 
and the “weights” A\, A 2 , . . . , A n . 

One major source of the previous fomula is in the theory of polynomial interpo- 
lation. If the nodes have been fixed, then there is a corresponding Lagrange interpo- 
lation fomula: 

P(x) = £/(*<)*••(*) 


1 = 1 


where L(x) = f] (— — — ). 




X\ X j 


This formula provides a polynomial p of degree < n — 1 that interpolate / at the 

nodes; that is, p(x,) = /(x t ) for 1 < i < n. If the circumstances are favorable, p 

will be a good approximation to /, and / f(x)dx will be a good approximation to 

J a 


j 


p(x)dx. Therefore 

f f(x)dxK [ p(x)dx = Y, /(x«) / h(*)dx = ^Aif(xi). 

Ja Ja Ja ■ i 


i=l 

rb 


t=l 


where we have put A, = f /,(x)<fx. From the way in which the previous formula 

J a 

has been derived, we know that it will give correct values for the integral of every 
polynomial of degree < n — 1. 

In the preceding discussion, the nodes have been arbitrary, although for practical 
reasons they should belong to the interval in which the integration is to be carried 
out. Gauss discovered that by a special placement of the nodes, the accuracy of the 
numerical intgration process could be greatly increased. There exist tables in which 
the numerical values of the nodes and weights are listed for the formula 


Those nodes Xj, X 2 , • • • , x n are called the Gaussian nodes with which the previous 
formula w'ill be exact for all polynomials of degree < 2n — 1. 

Here the table for n = 10 is given. 



x (n) 

x k 



± 0.14887 43390 

0.29552 42247 


± 0.43389 53941 

0.26926 67193 

n=10 

± 0.67940 95683 

0.21908 63625 


± 0.86506 33667 

0.14945 13492 


± 0.97390 65285 

0.06667 13443 


Usually, the integral interval of function /(x) is on [a,b]. In order to use the 
Gaussian quadrature formula, one should transfer interval [a,b] to interval [-1,1] using 
a linear transformation, for example, for the twisted-pair case in section 6.1, the 
following linear transformation can be carried out: 

/ f{0)dO = 7T / f(irt + (2 i -f 1 )r)dt 

J2vi J - 1 

where, 

0 nt -J* (2i -j- 1 ) 7r 
if 6 = 27rz , then t = — 1 
if 6 = 2t r(i +1), then t = 1 
d0 = 7 rdt. 

Now, one can use the Gaussian quadrature formula as follows: 

/(«)<» = »E + (2* + 1)') 

• /2,r ‘ fc=i 

For further details reference [12] should be consulted. 
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Chapter 1 


Introduction 


The goal of this research work is to analyze the potential effects of electromag- 
netic interference (EMI) originating from power system processing and transmission 
components for NASA Space Station Freedom. Particular emphasis is given to EMI 
sources radiating in a plasma environment, i.e., the ionosphere. The approach con- 
sists of three sequential steps: 

1. Determine the characteristics of the ionospheric plasma and formulate a work- 
ing model for purposes of EMI analysis. 

2. Develop analytical tools (models and computer programs) to predict the elec- 
tromagnetic fields from isolated sources radiating in unbounded free space or 
a plasma environment. 

3. Develop numerical tools (moment-method formulation) to predict the elec- 
tromagnetic fields from radiating sources on a simple space station structure 
in free space or a plasma environment. 
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1.1 Ionospheric Modeling 


There is considerable information on the properties of the ionosphere and the 
different natural phenomena that affect its electromagnetic characteristics. To start 
with, the ionosphere is a region of the earth’s atmosphere composed of molecules 
that are ionized by solar radiation. Figure 1.1 shows the molecular content of 



Figure 1.1: Atmospheric molecular density versus altitude [1]. 

the atmosphere as a function of altitude. Due to gravity, the heaviest molecules 
( 7 V 2 , O 2 ) are at the lowest altitude and the lighter molecules (He, H) are at higher 
altitudes. In the middle region, where the space station will be located, atomic 
oxygen (0) is the most prevalent. 

The degree of ionization of these molecules depends on many different factors. 
The two strongest contributors are time of day and altitude. During the nighttime, 
solar radiation does not reach the ionosphere and ions recombine resulting in the 
electron density, N ) decreasing. Figure 1.2 shows this effect and also the altitude 
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Figure 1.2: Time of day variations of electron density versus altitude [2]. 

variation. The envelope can be qualitatively explained in simple terms. At higher 
altitudes the solar radiation is very intense but there are few molecules. As the 
altitude decreases, the molecular concentration increases and more solar radiation 
is absorbed. Eventually, a point is reached where the high molecular concentration 
absorbs most of the remaining solar energy and the ionization reaches a maximum. 
Continuing to lower altitudes the electron density decreases. Since the atmosphere 
is composed of many different gases with different ionization and recombination 
characteristics, different layers and peaks in the electron density form [2]. These 
have been labeled the D, E, and F layers. An important consequence of the layering 
is that radio waves in certain frequency ranges can be bounced off of the ionosphere 
to improve communication distances. Figure 1.3 shows how some of the properties 
of the ionosphere are exploited in today’s world for useful means. 

Another important contributor to the electron density is sunspot activity. Fig- 
ure 1.4 shows sunspot activity observed for the past few cycles, and the predicted 
activity using the Lincoln-McNish method [3]. Increased sunspot activity increases 
solar radiation which leads to increased ionization of the atmosphere. This of- 
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Figure 1.3: Simple examples of ionospheric impact on radiowave systems [3]. 
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Figure 1.4: Observed and predicted sunspot numbers for several solar cycles [3]. 
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ten disrupts the normal communication schemes in positive and negative ways. 
On one hand, the layers become more pronounced which increases reception and 
transmission distances in radio communications. On the other hand, satellite com- 
munications can deteriorate if the frequencies are not high enough. 

At the space-station altitude range of 450 to 550 km, the ionosphere is considered 
fully ionized. It is a neutral region of mobile positive ions (predominantly 0 + ) 
and electrons (e - ) that form a plasma. Plasmas are considered the fourth and 
most energetic state of matter. They are also the most abundant state in the 
universe [1], and occur in everyday life in fluorescent lighting, vacuum tubes, video 
display terminals, and lasers. 

A unique feature of the ionospheric plasma is the presence of the earth’s mag- 
netic field. This makes the ionosphere an anisotropic plasma, where the polarization 
of the medium, P, and the electric field, E, are not in the same direction [4]. The 
Lorentz formula, F = g(E + vx B), states that a free electron traveling with a 
velocity, v, perpendicular to a constant magnetic field, B, will travel in a circular 
orbit. This leads to an angular (cyclotron) frequency for an electron in the earth’s 
magnetic field of 5 x 10 -5 Wb/m equal to 1.4 MHz. A consequence of this is that 
radiowaves close to 1.4 MHz are highly attenuated [2]. Also, the current density, 
J, and the electric field, E, are related by a tensor (matrix) quantity. It will become 
evident from the formulations in the next chapter that a tensor permittivity instead 
of a scalar value severely complicates the analysis of wave propagation in a plasma 
medium. 

In order to solve EMI type problems, the ionospheric plasma is usually dis- 
cussed in simpler terms. The zero-order plasma model is used, where the effects 
of the magnetic field of the earth are neglected. The theoretical formulation of 
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electromagnetic wave propagation is then based on a scalar permittivity. With this 
approximation, electromagnetic wave propagation in a plasma becomes similar to 
wave propagation in a waveguide. The zero-order plasma model also assumes a 
homogeneous plasma. Although the ionosphere is locally homogeneous, motion of 
the space station through the plasma creates ram/wake effects [5, 6] similar to that 
of an aircraft passing through air. Thus the ion and electron densities are increased 
in front of the space station and decreased behind the space station. 


1.2 Electromagnetic Analysis 

Computer analysis of EMI is a widely-used technique and sophisticated com- 
puter programs are available to simulate the EMI for complex systems involving 
hundreds of sources and susceptors [7]. Built into these programs are different types 
of models for sources and susceptors and transfer functions for various forms of cou- 
pling. To use the programs for analysis, a system must first be described in terms 
of sources and susceptors and their geometrical arrangement. This is then entered 
as parametric data into computerized data sheets. This procedure simplifies EMI 
analysis on major aircraft and weapons systems, where the programs have been de- 
veloped and proven effective. However, the complex procedure is not appropriate 
for modeling electromagnetic radiation in the ionospheric plasma. For this reason, 
attention turned to finding a general numerical code that could be easily worked 
with, understood, and modified to predict EMI in a plasma environment. 

The moment method is a numerical technique for solving electromagnetic field 
problems. It was first unified into a general procedure by Roger F. Harrington in 
1968 [8]. Since then it has been utilized on powerful computers to solve complex 
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electromagnetic radiation and scattering problems. The moment method utilizes 
a few assumptions that increases its ability to handle a problem as large as the 
space station. The primary assumption is that all current flows in thin wires and 
is axially directed. Both radiating sources and the space station structure are then 
modeled with wires, like the airplane in Figure 1.5. The result is that the radiated 



Figure 1.5: Wire-grid model of F-16 aircraft [9]. 


electromagnetic fields can be calculated anywhere in space relative to the structure. 
This shows what effect variations in the amplitude and location of sources on the 
space-station structure have on potential EMI problems. It should be pointed 
out that a parallel effort is under way to use the finite-element method for specific 
source modeling [10, 11] that will lead to the development of equivalent sources that 
can be input into the moment-method analysis. But, at this time, tools are being 
developed that are more general in their source modeling to allow for flexibility in 
design parameters. 

MININEC3 was the selected moment-method code [12] to be modified for two 
reasons. The first is that it is proven and well documented. It runs on a personal 
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computer, but it is limited by the size of the problem it can handle and by the speed 
at which it runs. This makes it ideal for the analysis of simple models, but not 
practical for more complicated and defined structures. The second reason is that 
MININEC3 is a scaled-down version of NEC, numerical electromagnetic code [13], 
which is one of the most advanced computer codes available. If MININEC3 can be 
modified for a plasma environment to accurately predict EMI for small problems, 
then the same methodology may be applied to NEC (or other moment-method 
codes) to handle more complicated problems. 


1.3 Overview of the Thesis 

Many steps are required before the electric fields in the ionospheric plasma 
medium can be predicted. Chapter 2 starts the work with theoretical formulations. 
First, the analytical solution for a radiating source is derived. Then the moment- 
method numerical technique is explained. Lastly, the zero-order plasma model is 
discussed. In chapter 3 analytical results are compared for the radiating source 
in free-space and a zero-order plasma. Chapter 4 starts with verification of the 
moment-method code in free space and zero-order plasma. Then the simple space 
station structure is explained and electric fields in free space and a zero-order plasma 
are presented and compared. Chapter 5 contains the conclusions which were based 
on the research and gives some suggestions for further work in this area. 
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Chapter 2 


Theory 


This chapter has three goals. The first is to derive the analytical solution for 
the electric fields from a radiating elemental dipole. The elemental dipole was cho- 
sen because its analytical solution is well understood and the theory behind the 
moment method is closely related to that of the elemental dipole. The second goal 
is to formulate and explain the moment method. The moment method is a numer- 
ical approximation for determining electromagnetic fields from complex thin-wire 
geometries. The final goal is to explain the zero-order plasma approximation. This 
includes a consideration of electromagnetic wave propagation and incorporation of 
the plasma into the moment method. 


2.1 Radiating Sources — Analytical Solution 


The formulation for determining the electric fields generated from an elemental 
(short) dipole are documented in many textbooks [4, 14, 15]. The following for- 
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mulation serves two important purposes. First, the analytical solution is used to 
formulate the moment method. Second, a closed-form solution is developed for the 
dipole radiating in a zero-order plasma. Since the elemental dipole is the fundamen- 
tal building block of this research, a brief description of it precedes the formulation 
of its fields. 

An elemental dipole is a short linear conductor. By definition, its length, L , is 
much less than the wavelength, A (L <C A), and the radius, a (d = 2a), is small 
compared to the length (a <C L ). The current, I, can then be assumed constant 
over its length. Figure 2.1(a) shows an elemental dipole and Figure 2.1(6) represents 
its equivalent. The transmission lines are balanced and shielded to provided the 



(a) ( b ) 

Figure 2.1: An elemental dipole (a) and its equivalent (6) [14]. 

uniform current and no radiation. They are therefore, excluded in the equivalent 
model. The elemental dipole can then be described as a thin linear conductor with 
a constant current. The electric field components E r , Eg, and E^, at an observation 
point, P, from a dipole coincident with the z-axis and at the origin, are defined in 
a spherical coordinate system as shown in Figure 2.2. 

The formulation for the electric field components begins with Maxwell’s equa- 
tions relating the electric field intensity, E, to the magnetic field intensity, H, 

V x E = —jujfxH., (2-1) 
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Figure 2.2: Electric field radiation of elemental dipole in a spherical coordinate 
system [14]. 

V x H = jue E + J, (2.2) 

where J is the source current density, uj is the angular frequency, /x is the perme- 
ability of the medium, and e is the permittivity of the medium. Using the identity 
V- V x (■) = 0 on (2.1), 

V • V x E = —jijJfi (V ■ H) = 0. (2.3) 

Since a vector with zero divergence can be expressed as the curl of some other 
vector, H can be written as the curl of A, 

H = V x A, (2.4) 

where A is defined as the magnetic vector potential, MVP. Substituting (2.4) into 
(2.1) gives 

V x E = — ju>/xV x A, (2-5) 

VxE + Vx ju>fi A = V x (E 4- ju>fi A) = 0. (2.6) 
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Since any curl-free vector is the gradient of some scalar, 

E + jufiA = — VV, (2-7) 

where V is the electric scalar potential. Substituting (2.4) and (2.7) into (2.2) yields 

V x V x A = — jweVV -I- k 2 A + J, (2-8) 

where k is the wave number, u>y/]Ie. Using the vector identity V x V x A = 
V (V > A) — V 2 A, (2.8) reduces to 

V(V- A) - V 2 A- Jfc 2 A = -jueVV + J. (2.9) 

Since only the curl of the MVP has been specified, its divergence can be chosen 
to simplify (2.9). Common choices are the Coulomb Gauge and the Lorentz Gauge. 
In this case the Lorentz Gauge 

V • A = -jojeV (2.10) 

is used to eliminate the electric potential term. Substituting this into (2.9) reduces 
the formulation to 

V 2 A + k 2 A = — J, (2.11) 

which is the Helmholtz equation or the complex wave equation. It relates the MVP 
to the source current density. Substituting (2.10) into (2.7) relates the electric field 
to the MVP. The result is 

E = —ju>fiA + t— V (V • A) . (2.12) 

jue 

The electric field resulting from a source current density can be determined by 
simultaneously solving (2.11) and (2.12). But this can be simplified further for an 
elemental dipole, by solving for the MVP and combining the two equations. As 
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previously stated, the elemental dipole is centered at the origin, and is coincident 
with the z-axis as shown in Figure 2.2. Since the current on the wire is in the z 
direction, the source current density and MVP are also in the z direction. At points 
away from the origin (not on the wire), J = 0 and (2.11) reduces to 

V 2 A, + k 2 A z = 0, (2.13) 


where A z denotes the (scalar) MVP in the z direction. For observation points 
sufficiently far enough away from the dipole, the source looks like a point. The 
resulting A z from a point source, Idz', is spherically symmetric (a function of the 
radius only), A z = A z (r). Evaluating the Laplacian of the MVP with the radial 
dependence only and substituting into (2.13) gives 


l_d_ 
r 2 dr 


, dA z 

dr , 


-f - k 2 A z — 0 


(2.14) 


in spherical coordinates, where r is the distance between the source and observation 
points. The two solutions for the MVP from (2.14) with r / 0 (not at the origin) 


are 


and 


a. = 


A, = — e jfcr . 


(2.15) 


(2.16) 


Since (2.15) represents an outward traveling wave and (2.16) represents an inward 
traveling wave, (2.15) is the appropriate solution for the radiating dipole. In order 
to find the constant C\, let u> — ► 0. Then (2.15) becomes 


A.-* 

T 


(2.17) 


and Helmholtz’s equation (2.11) becomes 
V 2 A, = -J z . 
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(2.18) 



Helmholtz’s equation for the case u> = 0 is solved by noting its resemblance to 
Poisson’s equation, 


V 2 V = , 

e 


where p v is the charge density. The solution of (2.19) is known to be 


( 0 = / 

J Ui 


pjjw 


voi 47 rer 


(2.19) 


( 2 . 20 ) 


where r', the source point, is the location of p v and r, the field point, is where V is 
observed. The volume of the source region is vol, and r = |r — r'| is the distance 
between the source and field points. Noting the mathematical analogy between 
(2.18) and (2.19), A z is determined to be 
Jz (r ; ) dt/ 


L (r) = / 

Jvc 


( 2 . 21 ) 


)voi 47r r 

Since the volume of the dipole can be written in terms of its cross sectional area 
times its length, L , and the cross sectional area times the source current density, 
J*, equals the source current, 7, (2.21) becomes 
Idz' 


[z ~ Jl 


'l 47t r 

Equating (2.22) and (2.17), results in 

Idz' 


C '=L 


L 47T 


Substituting back into (2.15) leads to the solution 


-k 


Idz’ 


-i kr 


L 47r r 


Performing the integration for the elemental dipole, 
A, = 

47rr 


( 2 . 22 ) 


(2.23) 


(2.24) 


(2.25) 


This shows that the MVP can be determined at a specified distance from an ele- 
mental dipole current source. If the magnetic field intensity were desired, it would 
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be determined by taking the del cross product (or curl) of the MVP, (2.4). This is 
the point from which the moment-method formulation will be started in the next 
section. 


Since it is common for the analytical solution of E from an elemental dipole to 
be expressed in spherical coordinates, the MVP in the z direction is transformed 
into spherical coordinates using 

A r = zA z • r = zA z • z cos 8 — A z cos 6 (2.26) 

and 

Ae = zA z • 9 = zA z • z (— sin 9) = —A z sin 8, (2.27) 

where A r is the radial component of the MVP, and A$ is the theta component of 
the MVP. The electric field expressed by (2.12) is now 

E = -juft (rA r + 9A e ) + -^V (V • (rA r + 8A e )) . (2.28) 

Evaluating V (V • A) in (2.28) and reducing E to the radial, E r , and theta, Eg, 
components yields 

1 d ( 1 d , . 1 d 


E r = -jufiAr + — ^ [ ^i~r 2 Ar + -^—z^A e sin9 


juie dr \r 2 dr 


r sm 


8d9 


and 


„ 4 1 d ( 1 d 2 . 1 d . . a 

E e = -]u>(iA e + - — ( — — r A T + — r- 7 — A* sm 8 


) 


(2.29) 


(2.30) 


jwer d8 \r 2 dr ~ r rsin 8 dd 

The phi component, shown in Figure 2.2 is not present. Substituting (2.26) 

and (2.27) into (2.29) and (2.30), and simplifying the expressions yields 

_ . , . cos 9 d ( 1 d / 2 a \\ 2 cos 6 d ( A z \ . . 

E, = -,u V A. cos * + — Fr (r A.) j - (— ) (2.31) 

and 


sin# d / 9 . \ . 2 A z sin# 


„ , . * sin (7 o ( 2 . \ 

E e = ju>fj,A z sm 9 — - — (r AA + — 

juter* or v ' 


juer A 


(2.32) 
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Substituting (2.25) into (2.31) and (2.32) and simplifying the expressions results in 
the final analytical solutions for the electric field intensities, 

IL 


E r = 


,-i* 


27T 


V e juter* ) 


cos 6 


and 


_ IL_ jkr Ijuw 


Eg = — e 
47 r 


+ +t-~] sin0. 

r V e ju>er* J 


(2.33) 


(2.34) 


Figure 2.3(a) shows Eg and Figure 2.3(6) shows E r . In three dimensions, Eg has 



Figure 2.3: (a) Eg and (6) E r for an elemental dipole source on the z-axis [14]. 

a characteristic donut shape around the dipole and E r appears as equally sized 
spheres on the top and bottom of the dipole. 

Equations (2.33) and (2.34) are often referred to as the near-field patterns for 
radiation from an elemental dipole. This is because, near the dipole E r and Eg are 
similar in magnitude. But far away from the dipole Eg E r , so that E r can be 
neglected. How far away is determined by examining the 1/r, 1/r 2 , and 1/r 3 terms 
in (2.33) and (2.34). The three terms are equal at r = A/2 tt in free space, where A 
is the wavelength. For r the 1/r term is dominant and the electric field 

reduces to 

Eg = ^ e - jkr ^ sin0. (2.35) 

4w r 
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If a frequency of 1 MHz (A ~ 300 m) is used as an example, the far field approx- 
imation could be used if r 50 m. Since the maximum dimension of the space 
station is about 100 m, the far field approximation cannot be used at this frequency. 
The attention now turns to formulating the moment method. 


2.2 Radiating Sources — Moment Method 

The moment method is a numerical tool for predicting fields from thin-wire 
structures for which there are no analytical solutions. It is a numerical method 
capable of connecting a large number of small wire segments and determining field 
patterns by enforcing electromagnetic boundary conditions on and between the seg- 
ments. Since it is based on current segments and the elemental dipole approximates 
a current segment, the analytical sol u tion developed in the previous section for the 
elemental dipole is utilized. This explanation of the moment method is similar to 
that in Kraus [14]. 

In the previous section, the MVP was related to the electric field E in (2.12). 
For an elemental dipole along the z-axis, it was shown that the MVP is in the z 
direction only. Substituting A z into (2.12) and evaluating V (V • zA z ) yields 
1 O' 1 A \ 

E ‘ = ]Zl “ iwftA - (2 ' 36) 

which reduces to 

E -=Mi£ +k ' A ‘)- (2 37) 

From the solution of the Helmholtz equation, the MVP was related to a current 
pulse segment in (2.24). Substituting the value for A z into (2.37) yields Pockling- 
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ton’s equation [14], 


E -M = ^J,{&’{^) + T'-’ kr ) I ^ d *' 


(2.38) 


where E z (r) is the radiated electric field at an observation point r from a current 
source at z' and r is the distance between the two points, r = |r — r'|. 


In the moment method, the solution in (2.38) is used twice. Obviously, the for- 
mulation may be used to calculate the radiated electric fields from current sources. 
But this is not chronologically the first time that it is used in the moment-method 
solution. First, (2.38) is used to obtain the current in the wire from the voltage 
applied across the dipole terminals. Then the radiated fields can be calculated. To 
explain this, a closer look at the dipole source is needed. 

The dipole source is physically a cylindrical current-carrying perfect conductor 
as shown in Figure 2.4(a). Assuming that the conductivity of the wire is high, the 
current is located entirely on the surface of the cylinder. For the dipole centered 
at the z-axis, all of the current is then at a wire radius a from the z-axis. This is 
shown in Figure 2.4(a). Since all of the current is equidistant from the center of the 
cylinder (z-axis), the current density will be considered to flow in a thin filament a 
distance a from and parallel to the z-axis. This is shown in Figure 2.4(6). 


In any electromagnetic radiation or scattering problem involving a conductor, 
there are two components of electric field, E app , which creates the current on the 
conductor, and E ra< j, which is the electric field radiated by the current carrying 
conductor. In the case of the elemental dipole antenna, E app results from a voltage 
applied across the antenna terminals (see Figure 2.1(a)). The total electric field 

E app -f E rad = E 0 b s , (2.39) 

where E 0 b s is the total observed electric field intensity. Since the skin depth is zero 
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Figure 2.4: (a) Cylindrical conductor model of elemental dipole. ( b ) Current fila- 
ment model used to replace the cylindrical conductor [14]. 


for a perfect conductor, the observed internal electric field, E 0 b s , is zero and 

E ra d = — E app . (2.40) 

For the special case of the radiated electric field within the elemental dipole, (2.38) 


becomes 


E a pp (z) — 


1 


d 2 




47 rjutc Ji \dz 2 
1/2 


k 2 

+ — e 
r 


-’’"'j I(z')dz’, 


(2.41) 


where r — [(z — z') 2 + a 2 j as shown in Figure 2.5. Since E app is known, (2.41) is 
an integral equation which can be solved for I (z 7 ). Then J(z 7 ) can be substituted 
into (2.38) to obtain the electric field everywhere. 


The electric field formulation is made more concise by substituting the Green’s 
function notation, 


G„. = 


,-jkr 


(2.42) 
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Figure 2.5: Current filament Idz' with field E t ( z ) at a distance r. 
into (2.41). The result is 

- JTJZ-J, {B a “' + k2 °“') < 2 ' 43 > 

It is shown in the Appendix that evaluation of the partial derivatives in (2.43) gives 
-E app (z) = + 3 kr ) ^2-3 j + k 2 a 2 ^j I(z') dz'.(2.44) 

Substituting u> = -^= and the intrinsic impedance of the medium r\ = yields 

-SwM = ^U~S~ ((! +^ r ) ( 2 - 3 (;)’) + lV ) ( 2 - 45 > 

The relationships in (2.38) and (2.45) form the theory behind the moment 
method. They are the relationships between applied electric fields, radiated elec- 
tric fields, and current sources. They have come directly from Maxwell’s equations 
with no assumptions about the medium that the radiation is traveling in. They are 
therefore, valid for free space or a plasma, using the proper permittivity and per- 
meability. The remaining part of this section incorporates the methodology behind 
the evaluation of (2.45). 
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The first transformation is to let observation points be denoted by subscript m 
and source points by subscript n. This will be done gradually over several steps 
starting with (2.45). It can be written as 

-E.„(,z) = J t nz')G(r mn )dz', (2.46) 

where 

G (r "" ) = wr ? ( (1 + ikT) ( 2 - 3 (;) 2 ) + ) - < 2 47 > 

and r mn is substituted for r, where r mn is the distance between source point n 
and field (observation) point m. The function, G(r mn ), consists of terms forming 
a geometry matrix of source and observation points that are used in the moment 
method. In any problem there are two of these matrices formed. One for the applied 
field to calculate currents in filaments, and another to calculate radiated electric 
fields from the currents. Although G (r mn ) may seem like a complicated term, in 
the moment-method formulation it is often simple compared to the current term. 
The environmental parameters are also contained in G(r mn ), but they do not vary 
from source to observation point. 

The current term in (2.46) can be approximated by the current series 

= (248) 

n=l 

where I n is a constant magnitude and F n ( z ') is a pulse function (equal to zero or 
unity) for incremental segments A z' n . In other words, the current is constant over 
a segment n in the same manner as with the elemental dipole. There are more 
complicated current series using overlapping segments with piecewise sinusoidal or 
triangular current distributions [14], but in MININEC3 the current pulse is used. 

The contribution from each current pulse segment n can then be summed to 
obtain the electric field at m. Incorporating this into (2.46) for the m th segment 
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results in 


- E app {z m ) = f; f I n F n {z')G{r mn )dz\ (2.49) 

n=l 

where the integral and the summation have been interchanged and A z' n is the 
length of the n th segment. The current portion can be pulled outside of the integral 
because its magnitude does not vary over a given segment. The result is 


-E app (z m ) = £)/„/ G(r mn )dz'. 

n=l JA < 


(2.50) 


Letting 


G mn = f G(r mn )dz’, (2.51) 

and assuming G(r mn ) is constant over the small interval A z„, (2.51) becomes 

Gmn = G (r mn ) Az' n . (2.52) 


(2.50) can now be reduced to a matrix form. 

- [E m ] = [G mn ] [/„] , (2.53) 

where [2? m ] is the matrix of electric fields from the matrix of current segments, [/„]. 
Multiplying both sides by the segment length A z, 

~ [V m ) = [Zmn] [/n] , (2-54) 

where V m is the applied voltage at the m th segment, and Z mn is the impedance 
matrix. This equation can be solved for the vector of current amplitudes [/„]. These 
can then be substituted into (2.38) to obtain the radiated electric field anywhere. 


2.3 Ionosphere — Zero-Order Plasma 


The ionosphere is approximated by neglecting the earth’s magnetic field. This 
is the zero-order plasma approximation. In this section, electromagnetic waves 
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radiating from an elemental dipole in a plasma will be shown to behave like an 
electromagnetic wave propagating in a waveguide. The formulation parallels the 
explanation given by Collin [2], At the end, the zero-order plasma model is incor- 
porated into the moment-method formulation. 


It seems appropriate to start the formulation for the zero-order plasma in simple 
terms. Considering a moving electron point charge in the presence of an electric 
field, E, the force on the particle is given as 

m e ^ = -eE, (2.55) 

where m e is the mass of the electron, v its velocity, and — e its charge. For the 
sinusoidal case, (2.55) becomes 

jcum e v = — eE. (2.56) 


Since the current density, J, equals eNv, substituting this into (2.56) yields 

Ne 2 

J = —eNv = - E. (2.57) 

ju>m e 

This current density can then be substituted into Maxwell’s equation (2.1), 

Ne 2 


V x H = ju;e 0 E + J = ju)e 0 I 1 


ujm. 


E, 


(2.58) 


where e 0 is the permittivity of free space, 8.854 x 10 12 (F/m). The equivalent 


permittivity of the zero-order plasma is then 

e = e 0 (l ~— V (2.59) 

\ u>m e J 

which is the scalar quantity times e 0 mentioned in chapter 1. 


Plasma models can be discussed in more useful terms using their characteristic 
plasma frequency. It is defined as 


w p = 


\N£_ 

m e e 0 ' 


(2.60) 
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Writing e in terms of u p yields 



Since the wave number, k , is defined as Wy/eji, k becomes 

k = 1 - (^f) • 


(2.61) 


(2.62) 


The permeability of free space, is introduced since the plasma is not a magnetic 
material. Writing k in terms of wavelength, A, 

27T 


k = 


f l 


-(?)' 


(2.63) 


In section 2.1, the analytical solutions for the electric field components from the 
elemental dipole were shown to be 


_ IL •. ( \ [a 1 \ 

E r = — e 3 —J~ + “ 3 cos ^> 

Z7T \r* V e juier* ) 


and 


Eg = !±e~”" (t£ 


+ + ■ 1 3 ) 
r i V e juer* J 


(2.64) 


(2.65) 


47t \ r 

in (2.33) and (2.34). Substituting the expression for k from (2.63) into the expo- 
nential term yields 


e-*' = e~’ rlf ' 


( 2 . 66 ) 


In the zero-order plasma approximation, k is either a real or an imaginary term 
depending on the frequency ratio. Then for k = /3 — j a, where a equals the 
attenuation constant and (3 equals the propagation constant, either f3 = 0 or a = 0. 
Therefore, for u> > w p , the wave propagates, and for u> < u) p , the wave does not 
propagate, but decays exponentially as an evanescent wave, 


(2.67) 
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where 



( 2 . 68 ) 


This behavior in the zero-order plasma is similar to electromagnetic waves in a 
waveguide, where waves below a cutoff frequency become evanescent. 


It will be shown in the following chapter that the plasma frequency of the iono- 
sphere is typically in the 2 to 10 MHz range. Since the frequencies and harmonics 
from the space station power system processing and transmission components are 
below this, the radiated electromagnetic waves are evanescent. Electric field inten- 
sities in a zero-order plasma can then be determined by multiplying the electric 
field intensities in free space by an attenuation factor, e” ar , where a is the real 
number obtained from (2.66) when u) < u> p . The free space phase term can sim- 
ply be neglected since it does not effect the magnitude of the intensity. For the 
moment-method formulation, e~ ar can replace the e“ jfcr term in (2.38). 

It is important to summarize the approximations that have been made in ac- 
counting for the plasma: 


1. the earth’s magnetic field has been neglected, therefore the ionospheric plasma 
is assumed to be isotropic (the zero-order plasma approximation); 

2. the frequencies of interest are assumed to be below the (electron) plasma 
frequency, hence evanescent rather than propagating fields result; 

3. the attenuation of the near-field evanescent fields is assumed to occur in the 
same manner as for the far fields; 

4. the attenuation factor, e” ar , is incorporated into the moment-method solution 
only in the calculation of the radiated fields and not in the calculation of the 
currents which produce the fields. 
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The approximations listed here simplified the development of a basic framework 
for EMI analysis using the moment method in the ionospheric plasma. Continuing 
work will improve the sophistication of the model by removing the first, third, 
and fourth approximations, thereby resulting in a more realistic simulation of the 
radiation problem. 
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Chapter 3 


Results — Analytical Solution 

In Chapter 2, the analytical formulations for E r and Eq from an elemental dipole 
were given in (2.33) and (2.34). In the last section of the chapter, the zero-order 
plasma model was incorporated into the formulations. In this chapter, the electric 
fields are calculated and then compared for unbounded free space and the zero- 
order plasma. But first, some input conditions are discussed as they relate to the 
space station and the ionosphere. Then computerized simulations are performed 
for different circumstances. 

3.1 Attenuation Factor 

The distinguishing property of the ionosphere was determined in Chapter 1 
to be the electron density. In order to determine the electron density, curves for 
the maximum and minimum charge density N, #/m 3 , as functions of altitude are 
shown in Figure 3.1 [5]. For a given altitude the maximum electron density occurs 
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Figure 3.1: Ionospheric charge density as a function of altitude [5]. 


during the daytime and the minimum occurs at nighttime (also see Figure 1.2). 
The two electron density curves in Figure 3.1 were digitized by approximating 
them with linear equations and entering them into the computer program. For a 
given altitude, a value is entered and used to interpolate between the minimum and 
maximum electron density values, i.e., unity corresponds to maximum N. In a more 
specific application, the program IRI-86 [6] could be used to give a more accurate 
value of N based on all the parameters that influence the ionosphere. Since the 
plasma serves to attenuate EMI, the worst case is for minimum N or free space. 
Meanwhile, it should be pointed out that overdesigning EMI protection for the free 
space case could be costly. 


The space station is being designed to operate between an altitude of 450 and 
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Figure 3.2: Plasma frequency versus electron density. 

550 km. A conservative range for N would then be between 5 X lO 10 and 2 x 10 12 
electrons per cubic meter [5]. Figure 3.2 shows the plasma frequency versus the 
electron density for the desired altitude range calculated from (2.60). Similar to 
a waveguide, propagation occurs for frequencies above the plasma frequency and 
attenuation below the plasma frequency. From the graph, it can be seen that the 
minimum plasma frequency of 2 MHz corresponds to the minimum electron density. 
As stated in Chapter 1, this analysis concentrates on frequencies below the plasma 
frequency. Figure 3.3 shows this by plotting the attenuation versus frequency (2.68) 
up to the plasma frequency for given electron densities. The attenuation in dBjm 
is calculated from 201og 10 e a . The graph also shows that the dB/m attenuation is 
fairly constant over frequency variations up to near the plasma frequency. 
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Figure 3.3: Attenuation versus frequency for specified electron densities 
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3.2 Computer Simulation 


The computer program was written to calculate the analytical solutions for the 
elemental dipole in unbounded free space and the zero-order plasma. A quick review 
of the equations and parameters follows: 


Er = + 


l 


27 r 


cos 6 


]u>e 0 r* 


(3.1) 


and 


E, = IL e -*' (m ± + sin 0 

47 r \ r r| e 0 ju>e 0 r* J 


(3.2) 


As discussed in section 2.3, the free space permittivity and permeability are used 
in the formulation. The remaining parameter, k, was determined to be equal to 
u)^JpL 0 e 0 for unbounded free space and for the zero-order plasma 

approximation. For the plasma, k is either real or imaginary, giving propagation or 
attenuation, depending on whether the frequency, w, is above or below the plasma 
frequency, u> p , respectively. 


By substituting the parameters into E r and Ee, the electric fields are obtained 
as functions of (u,N,r,8). The current, /, and length, /, were set equal to unity, 
and for free space the phase term is neglected since only the magnitude of the 
electric field is desired in the case of the elemental dipole. The final result of the 
simulation is the electric field intensities for an elemental dipole in free space and 
a zero-order plasma at given points. 
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3.3 Zero-Order Plasma Compared to Free Space 


The following simulations look at the distinguishing characteristics of the electric 
field in the zero-order plasma compared to free space. The first case examines the 
electric field components as the distance from the source increases. The next case 
looks at the effect of changing the frequency. The final case addresses the issue 
of variations in the electron density. Before any calculations can be made, the 
parameters that do not vary for a particular simulation need to be set. The first 
two are the length of and current through the elemental dipole. It was previously 
mentioned that these are set equal to unity to normalize the electric fields. In the 
next chapter when MININEC3 is compared to analytical solutions, the current and 
length will be adjusted to physical values. The next consideration is choosing a 
value for 9. Since the effect of 9 on field components is given in (3.1) and (3.2), its 
value was set to 45° to normalize its contribution to E r and E$. 

To examine the plasma’s effect on the electric field as the distance from the 
source increases, it seemed appropriate to start the calculations at 1 m, and extend 
out to the present design length of the space station, 100 m. Figure 3.4 shows the 
fields over this distance for free space and the zero-order plasma. The frequency was 
chosen to be 20 kHz because it corresponds to the switching frequency of converters 
on the solar panels. The initial electron density, 10 11 #/m 3 , was chosen because it 
is in the middle of the specified range. As expected, the electric field components in 
the plasma are attenuated. The electric field is expressed in dBfiV/m to conform to 
NASA’s EMI documentation requirements [16]. This is calculated by expressing E 
in fiV/m, and taking 201og 10 E. It is effectively a comparison of the field compared 
to a unit value of 1 fiV/m. 
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Figure 3.4: Comparison of electric field magnitude versus distance for free space 
and a zero-order plasma environment. 



To examine the plasma’s effect on the electric fields at different frequencies, a 
range must be set. Remembering that waves propagate above the plasma frequency 
and are evanescent below, it wouldn’t make any sense to choose a value above the 
plasma frequency. What is known is that only values in the kHz range have ever 
been mentioned by NASA for some of the power system components. With this 
in mind, 1 MHz is established as the maximum test point. For the lower limit, 
the 20 kHz switching frequency is chosen. Figure 3.5 compares Eg for the different 
frequencies in both free space and zero order plasma. For both mediums, as the 
frequency increases the magnitude decreases. This was also noticed in the E T 
component. As was predicted in Figure 3.3, the dB attenuation in the plasma is 
similar for both frequencies. 

The final issue looked at is to vary the electron density. From Figure 3.3, an 
increase in N should increase the attenuation. Figure 3.6 shows Eg for a typical 
range of free space, N = 0, up to a maximum N = 10 12 #/m 3 . As expected, an 
increase in N increases the attenuation. 


3.4 Summary 

In this chapter, a method was presented for determining the attenuation factor 
from measured charge density curves. The dB fm attenuation from the zero-order 
plasma was shown to be constant up to approximately one-half the plasma fre- 
quency. The electric field components for the elemental dipole in free space and the 
zero-order plasma were compared as functions of distance, frequency, and electron 
density. 
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Figure 3.5: Comparison of electric field intensity for different frequencies, FS=free 
space, P=plasma. 
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Figure 3.6: The effect of varying the electron density on the electric field. 
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Chapter 4 


Results — Numerical Solution 


In Chapter 3, the electric field for an elemental dipole was compared for free 
space and a zero-order plasma environment. The goal of which was to obtain 
preliminary results for some simple cases. However, the analytical approach is 
only applicable to isolated sources, i.e., sources that radiate in unbounded space 
without the presence of other structures. In this chapter, results for a more realistic 
space station configurations are discussed. But first, a verification procedure for 
MININEC3 is presented. This is done to insure that the program produces correct 
results in both mediums, and to gain experience with choosing proper convergence 
parameters. 

MININEC3 is by no means a complex implementation of the moment method 
like NEC is. It is small, runs on a PC, is written in BASIC, is not cluttered 
with built in functions or designs, is well documented, and is easy to work with. 
It provides a good test bed for ideas. As MININEC3 is a widely accepted code, 
the validation procedures used to make sure the program produces correct results 
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center around the type of models needed for the space station structure, instead of 
complex antenna designs. This means taking simple dipole and monopole sources 
on simple wire structures and at relatively low frequencies in the KHz range. This 
way the results can be qualitatively understood and provide proof of concept for 
modifying a larger program like NEC. 

MININEC3 is validated using two techniques. How accurate are MININEC3 
results compared to analytical solutions? For problems were there are no analytical 
solutions, how does MININEC3 results compare to more complex moment-method 
solutions? Four tests were used to answer these concerns. The first test compares 
solutions for the input impedance for some simple antenna configurations. The sec- 
ond test looks at currents on a long-wire structure for different frequencies. A third 
test looks at the numerical accuracy for determining fax fields in free space from 
dipoles of different numbers of segments. The final validation technique compares 
MININEC3 and the analytical solution for near field radiation in free space and 
the zero-order plasma. In all four cases, convergence of the solutions is obtained by 
paying strict attention to the number of wire segments that are used to discretize 
a structure. A general rule of thumb is 10 segments per wavelength. This allows 
for accurate representation of the current waveform on the wire. 


4.1 Input Impedance 

The input impedance is a particularly sensitive measure of the numerical ac- 
curacy of MININEC3 [14]. The reason for this is the use of pulse functions for 
the current instead of more complicated functions (explained in section 2.2). In 
the first two input impedance cases, MININEC3 is compared to solutions using 
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Table 4.1: Input impedance of an isolated dipole, in ohms. 


N, Modes 

MININEC 

King's Solution* 

1 

N/A 

2.000-jl921 

3 

1.927-J1980 

1.892-jl916 

5 

1.902-jl970 

1.864-jl905 

7 

1.887-jl957 

1.856-jl899 

15 

1.856-jl922 

N/A 


* King calculated the solutions using piecewise sinusoidal functions [14]. 
length = A/10, wire radius = A/10000 

piecewise sinusoidal functions. In the third test, the analytical solution is known. 
Table 4.1 shows the input impedance for an isolated dipole antenna of length A/10. 
For convergence, the results show that approximating the antenna by only three 
segments is sufficient to obtain the correct input impedance. Table 4.2 shows the 
antenna impedance when another identical wire is located parallel to and at a dis- 
tance of A/ 100 from the first antenna. The wire coupling configuration is more 
representative of the type of problem that would be of concern for an EMI analysis 
on the space station. Once again, only a few segments are required to adequately 
approximate the antenna input impedance. Table 4.3 completes the study of input 
impedance by calculating it for a longer antenna, a A/2 dipole. Again, only a few 
segments are required to accurately approximate the input impedance. The rule 
of thumb dictated 10 segments/A for accurate results. For this dipole, 5 segments 
would then be needed, so these results justify the rule of thumb. < 
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Table 4.2: Input impedance of dipole 1 in the presence of dipole 2, in ohms. 


N, Modes 

MININEC 

King*s Solution* 

1 

1.264-jl908 

1.382-jl822 

3 

1.281-j 1918 

N/A 

15 

1.256-j 1868 

N/A 


* King calculated the solution using piecewise sinusoidal functions [14]. 
d=dipole separation distance= A/100 
length = A/10, wire radius = A/10000 


Table 4.3: Input impedance of A/2 dipole, in ohms. 


N , Modes 

MININEC 

Known[ 14] 

3 

72.39-i-j60.75 


5 

75.90H-j41.88 

73.0+j42.5 

9 

78.05H-j41.88 


19 

79.28+j42.09 



length = A/2, wire radius = A/10000 
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4.2 Current Distribution 


Another method used in validating MININEC3 was solving for the induced 
currents in the wire structure. These currents form a special type of conducted 
EMI, i.e., currents induced in the space station structure by EMI sources. As 
a test case, Figure 4.1 shows a 100 m wire (divided into 34, 3 m pulse segments) 


z 



x 


Figure 4.1: Simple wire model of the space station with a short monopole antenna 
located at one end. 

excited by a short monopole antenna at one end. The current distribution along the 
simple structure is then calculated for a 20 KHz , 100 KHz, 1 MHz , and 5 MHz 
source and shown in Figure 4.2. At the lower frequencies of 20 and 100 KHz , 
the normalized current distribution is identical. This is attributed to the length of 
the structure being less then 3 % of the wavelength. As the frequency increases 
to 1 MHz , the distribution changes only slightly as the structure is then 1/3 of a 
wavelength. At 5 MHz, which is higher than the maximum frequencies discussed, 
the structure is 5/3 wavelengths. This higher frequency was used because it provides 
nice qualitative results since the current standing wave pattern shows a wire length 
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Figure 4.2: Current distribution along a 100 m wire divided into 34 pulse segments 
with a monopole source on one end at 20 KHz, 100 KHz, 1 MHz, and 5 MHz. 

of 5/3 wavelengths as expected. If the results for 10 MHz were shown, the current 
waveform would approximate 3- wavelengths. In the 4 cases shown, the current 
solutions converge. 


4.3 Comparison of Far-Field Results 


Although the far-field approximation is not used for the space station config- 
uration, it is the next logical step in validating MININEC3 and further increases 
confidence in the program. In Figure 4.3, far-field radiation patterns for a A/2 
dipole antenna with an increasing number of segments calculated with MININEC3 
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Figure 4.3: Calculated far-held radiation pattern for A/2 dipole antenna for an 
increasing number of segments using MININEC3 and analytical solution. 

are compared to the analytical solution. As with the input impedance of the A/2 
dipole, the results converge to known values for 5 segments. The results for this 
graph were normalized by dividing through by the source current and length of the 
dipole. 


4.4 Comparison of Near-Field Results 


The final step in validating MININEC3 is to compare analytical near field solu- 
tions to MININEC3 solutions. This point had been foreseen and was a reason for 


43 




developing the dipole simulation. Instead of comparing MININEC3 to analytical 
solutions for each environment separately, it makes more sense to present them 
together. In that respect, a detailed description of the dipole and conditions used 
to determine the electric fields follows. 

The dipole used in MININEC3 approximates the elemental dipole used in the 
analytical solution. The elemental dipole assumes a constant current across its 
entire length. In the moment method the current was shown to be distributed over 
several pulse functions. Numerically this has a big advantage because the boundary 
conditions at the ends of the dipole are satisfied by forcing the currents to be zero. 
This was shown in Figure 4.2 but not pointed out until now. The elemental dipole 
on the other hand, had a constant current over its entire length. This difference 
needs to be corrected for comparison. The dipole used in MININEC3 is divided 
into 5 pulses with the center being fed by a 1 volt, 100 kHz source. The current 
distribution on it appears as a step-like triangular function. In a qualitative sense, 
the analytical current distribution would look like a triangular step distribution if 
half its length with zero current was added to each end. Therefore, the length used 
in MININEC3 needs to be twice as long as its equivalent elemental dipole. In both 
cases, the dipole is centered along the z-axis and the electric fields are calculated 
along the x and z axis. The x and z components in this configuration then translate 
to E t and E$ as formulated in (2.33) and (2.34). 

Figures 4.4 and 4.5 show that MININEC3 results converge very closely to 
analytical solutions. The wire length and source current determined by MININEC3 
(shown on graph) are used as inputs for the analytical solutions. This was done 
to obtain the actual electric field intensities in dBfiV/m instead of normalizing the 
data to the source parameters. As discussed in Chapter 2, the zero-order plasma 
fields are calculated by multiplying the contribution from each pulse segment by the 
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Figure 4.4: x-component of electric field calculated analytically and numerically for 
unbounded free space and zero-order plasma. 



Figure 4.5: z-component of electric field calculated analytically and numerically for 
unbounded free space and zero-order plasma. 
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attenuation term, calculated from frequency, electron density, and distance between 
source and observation point. The same electron density used previously is used 
here, N = 10 11 #/m 3 . 


4.5 Electric Field Contour Plots for Space Station 
Structure 

A simple wire model of Space Station Freedom is shown in Figure 4.6. The 100 m 



Figure 4.6: Wire model of NASA Space Station Freedom. 


longitudinal wire consists of 20 segments, with a one-volt source on the left-most 
segment. The source location was selected to simulate the switching associated 
with the photovoltaic solar panels. The two vertical wires, representing habitation 
and laboratory modules, are 50 m long and spaced 5 m either side from the center 
of the 100 m wire. Each 50 m wire is divided into 10 segments. Using the rule 
of thumb, 10 segments per wavelength, the shortest A for the 5 m wire segments 
would be 50 m, or / < 6 MHz. 
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A contour plot of the electric field intensity for unbounded free space is shown 
in Figure 4.7. To form this, the following parameters were used in MININEC3. The 
above mentioned space station model was entered in the x-y plane with the source 
on the 100 m boom at the origin and extending out into the positive x-axis. The 
free space environment was selected. A source with a frequency of 1 MHz was set, 
because it offers the worst case for attenuation when the zero-order plasma model 
is used. Starting at x = — 10 m and y = —35 m, near electric fields were calculated 
in 2 m increments till x = 110 m and y = 35 m. There is then a near-field moment- 
method approximation since multiple field points exist along a single source point. 
The contour lines are specified in 5 dByVfm increments. 

The results show the field strength to be relatively strong around the source as 
expected. The strong fields at the far end of the structure indicates that conducted 
EMI (currents induced on the structure) should be a major concern to the space 
station designer. Especially if the attenuation over distance is a main strategy em- 
ployed to alleviate EMI. The field then decreases as the distance from the structure 
increases to a value of about 50 dBfiV/m at the outer boundary of the grid. For the 
zero-order plasma the attenuation at the boundary is expected to be much greater. 

The same conditions used in the previous free-space contour plot were used to 
determine the fields in the zero-order plasma model. The only difference being 
the zero-order plasma environment was selected with an electron density of N = 
10” #/m 3 . A contour plot of these fields is shown in Figure 4.8. 

Immediately, the increase in the number of contour lines shows the electric 
field to be much more attenuated. The yellow lines starting at 45 dB(iV/m and 
below, were not on the free space contour plot. This is the principal benefit of 
the zero-order plasma, i.e., moving a platform a relatively small distance from the 
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Figure 4.7: Electric field intensity in dB/iV/m from radiating source on space 
station model in a free space environment. , 
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Figure 4.8: Electric field intensity in dB^iV/m from radiating source on space 
station model in a zero-order plasma environment. 



space station provides a cleaner electromagnetic environment with the plasma then 
free space would. The fields close to the source, less than 10 m, are only slightly 
attenuated as the theory would have suggested. But some of other positions near 
the structure provide some interesting results that were not expected. 

The plasma appears to limit the conducted EMI problem, but in doing so the 
fields near the center of the structure are increased. In some areas (around x=40 m 
and y=0 m) the field increases by around 30 dByLVjm. At first an error was 
suspected in the data entry because stronger fields in the middle of the structure 
would be caused by higher currents in the structure. The current data did not 
support this conjecture because the current in each pulse segment was identical in 
both mediums, as expected. So how could the sum of the electric fields from each 
pulse increase in the zero-order plasma? Looking back at the formulations and 
the difference in the two media indicates a possible explanation. In the zero-order 
plasma, the fields are evanescent and have no phase term. This makes the individual 
fields from each pulse segment add as scalar quantities. In the free space case, the 
fields have phase terms associated with them so they add as vectors. Therefore, at 
certain locations the field strength in the zero-order plasma could be higher because 
of field components out of phase negating each other in free space. This effect was 
not anticipated because the modeling of the plasma by attenuation has only been 
applied to the far field in the past. 


4.6 Summary 

MININEC3 has been shown to produce correct results in free space and the zero- 
order plasma by comparing solutions for input impedance, current distribution on a 
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wire structure, far-field patterns, and near-field patterns to solutions obtained from 
analytical results or more advanced moment-method codes. Results are presented 
in the form of contour plots for electric fields from a simple space station structure, 
showing the attenuation from the zero-order plasma compared to free space. 
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Chapter 5 


Conclusions and Further 
Research 


The contributions of this research work fall into two general categories. The first 
is in modeling of the ionosphere and developing a working model for electromagnetic 
wave propagation. The two interests that distinguish this area from others is the 
indepth analysis of frequencies below the plasma frequency and the interest in the 
near-field components. The second contribution centers around incorporating the 
zero-order plasma into the moment-method and ultimately generating field plots 
for a simple space structure. 

A large portion of this research work, as with any research work, involved col- 
lecting data, formulations, and theories on electromagnetic wave propagation in the 
ionosphere. At the altitude of the space station, the ionosphere was approximated 
as a zero-order plasma, i.e., the earth’s magnetic field was neglected. The character- 
istic properties are principally determined by the electron density which is mainly 


52 


a function of altitude, time of day, and solar activity. From the electron density 
the plasma frequency was determined. For radiation above the plasma frequency, 
the waves propagate. For radiation below the plasma frequency, the waves decay 
exponentially, i.e., they axe evanescent. This behavior is similar to electromagnetic 
wave propagation in a waveguide. It was also shown that up to about half the 
plasma frequency, the attenuation from the plasma is fairly constant. For the space 
station, half the minimum plasma frequency is about 1 MHz . Since the power 
system components and their harmonics are usually well under 1 Afi7z, it was set 
as the upper limit of frequencies to be considered in this analysis. 

The moment method was then formulated as an analysis tool for performing 
calculations. To do this, analytical solutions for the elemental dipole were developed 
for both media. The solution for the fields radiated by an elemental dipole provided 
the basis for an explanation of near and fax fields. This raised two important 
conclusions for the research. The first was that the far-field approximation could 
not be used for the space station versus the more complicated near-field solution 
because of the low frequency of radiation being analyzed. The second was that the 
plasma models that have been developed have only been applied to the far field, 
since the near field was never an issue before. The analytical solutions were then 
used to test the modified version of MININEC3. This computer program, chosen 
mostly for its simplicity in implementing the moment method, was shown to be 
a good program to do test-case calculations. The dipole results showed that the 
attenuation in a plasma medium was as expected and it provided accurate results 
for both free space and the zero-order plasma. 

Wire models were then formed for the space station. The main reason such a 
simplistic program could handle something as large as the space station is the low 
frequency consideration. In order for the moment method solutions to converge, 
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pulse segments on the structure should be less then 1/10 of a wavelength. The 
space stations longest part is the main boom at about 100 m. At the maximum 
frequency of 1 MHz (wavelength=300 m), the structure could be modeled with 
under 4 segments. But because of resolution for determining grid plots, a smaller 
segment length of 5 m was used in the structure. This was determined to be 
adequate for the resolution in the contour plots and also near the maximum size of 
a problem MININEC3 is designed for. Although a larger problems could be done, 
the time to calculate field patterns starts to become an issue. The first test was to 
calculate the current on a 100 m wire for several different frequencies. 

A more complicated model of the space station was then developed for calcu- 
lating field patterns. The electric field contour plots from the wire model for the 
two media raised several interesting points. The first was the relatively high fields 
at the opposite end of the structure resulting from the conducted EMI or currents 
induced in the structure. Also, from the zero-order plasma contour plot, fields were 
shown to be strongly attenuated for only a short distance from the structure com- 
pared to the free space case. A not-so-obvious advantage of this is that EMI from 
the space station will be reduced for structures that may be near the station, i.e., 
satellites and the space shuttle. 

The zero-order plasma plot also showed some increases in field intensities near 
the center of the structure. One conclusion was that this effect could be attributed 
to the lack of a phase term in the waveguide model. But in interpreting these 
results, they can only be as good as the model. Since this model is only proven in 
far field cases, where phase would not have a large effect on the field intensities for 
free space calculations, the waveguide model in the near field needs more work. It 
is thought that since the currents were the same on the structure in the two media, 
the fields are not expected to increase in a plasma. But the arguments seem to be 
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all qualitative in this area. 


Indications axe that the moment method, MININEC3 in this case, properly 
models the space station with thin wires and predicts the correct field results in 
both environments. 


5.1 Future Work 

There are many areas of work that need further research. The main focus should 
be around developing better models for near-field electromagnetic wave propaga- 
tion in the plasma surrounding the space station structure, In particular, adding 
the contributions for the electric fields from each pulse segment as scalar quantities 
should be addressed. More testing and development on structures with MININEC3 
needs to be completed to determine the effect of plasmas. The major problem of 
work in this area will continue to be in developing models for near-field electromag- 
netic radiation and analytical solutions for testing purposes. 
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Appendix A 


Evaluating Partial Derivatives 
and Simplifying 


The applied electric field at location z in terms of a current source at z' was 
given in (2.43) as 

~E app {z) = I{z')dz\ (A.l) 

where 

p-jkr 

G zz . = 

r 

and 

r = [(z — z) 2 + a 2 ] / . 

Evaluating in terms of r 

d 2 d 2 e- jkr 

dz* Gtz ' ~ dz> r ‘ 
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Taking one partial 

a £*; -jk rj-* (|) - e-*_ (|) 

dz r r 2 


5 e~ jkr -e~ jkr (dr' 


dz r 


r 2 V dz . 


(jkr + 1) . 


Then taking the second 


= ;V [ikt-‘ k (t) 2 (jkr + 1 ) - (If) (Jkr + 1) - 

e ~’ lr (H'fjk] + (£)’ (jkr + 1) 2r. 


= *7r jkr 1 (jkr + 1) (|f) ! - r ! (jkr + 1) gf- 

i^(|) 1 + 2r(j* : r + l)(|) J l. 


il G = 

dz 2 ° zz ’ 


C ^[-k 1 r 3 {^y-jkr^-r^+ 
2 jkr 3 (U) I + 2r- (fj) 2 ! . 


[<|> ! + 2 ^ ! + *) - £ (i^ + ^ 


:G„. = 


,-jkr \ ( Q r ' 


( -fc 2 r 2 + 2jkr + 2 ) - ^ (j&r 2 + r 


Evaluating the partial of r, 


= v/( z “ z ') 2 + ° 2 > 


— = (z - z')J(z - z') 2 + a 2 , 
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2 



(f - fT 

(z - z') 2 + a 2 ' 


But since (z — z 7 ) 2 = r 2 — a 2 , 



The second partial of r is 


d 2 r 
dz 2 


((z - z') 2 + a 2 ) 3 4- ((z - z') 2 + a 2 ) 2 . 


d 2 r 

dz* 


—r 3 (r 2 — a 2 ) -fir 1 . 


d 2 r a 2 — r 2 1 a 2 

dz 2 r 3 r r 3 


(A.3) 


5 2 r 1 /a 2 \ 

dz 2 r 2 y r ) 


(A.4) 


Substituting (A.3) and (A.4) into (A. 2) results in 

§^2 Gzz ' = K 7 " 2 “ ° 2 ) ( _fc2r2 + ^ kr + 2 ) - ° 2 (i^ r + 1)] • 

^-zG zz < = 6 fa 2 — jkra 2 — k 2 r A + 2jkr 3 + 2r 2 + a 2 k 2 r 2 — 2 jkra 2 — 2a 2 ] . 
az 2 r 5 1 1 


B 2 r -1 

— — = — — [— fc 2 r 4 + 2jkr 3 + 2r 2 + a 2 k 2 r 2 — 3jfcra 2 — 3a 2 ] , (A. 5) 

dz 2 r b 1 J 

The final substitution is (A. 5) into (A.l). 


- E °pp W = sfc ft *3r (-*V + + 2r ! + 

—3 jkra 2 — 3a 2 + & 2 r 4 ) 7 (z 7 ) dz 7 . 

—E app (z) = — — — ■ — J — y~ (+2j fcr 3 + 2r 2 + a 2 k 2 r 2 — 3j fcra 2 — 3a 2 ) 7 (z 7 ) dz 7 . 

-E app (z) = 4 ^- we ^(1 + jkr ) ^2 — 3 ^ + fc 2 a 2 ^ 7 (z 7 ) dz'. (A. 6) 

(A. 6) is the simplified form referred to in the text. 
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